2#include <pinocchio/algorithm/frames.hpp>
3#include <pinocchio/algorithm/joint-configuration.hpp>
4#include "robot_model/Model.hpp"
5#include "robot_model/exceptions/FrameNotFoundException.hpp"
6#include "robot_model/exceptions/InverseKinematicsNotConvergingException.hpp"
7#include "robot_model/exceptions/InvalidJointStateSizeException.hpp"
10Model::Model(
const std::string& robot_name,
const std::string& urdf_path) :
11 robot_name_(std::make_shared<
state_representation::Parameter<std::string>>(
"robot_name", robot_name)),
12 urdf_path_(std::make_shared<
state_representation::Parameter<std::string>>(
"urdf_path", urdf_path)) {
17 robot_name_(model.robot_name_),
18 urdf_path_(model.urdf_path_) {
23 std::ofstream file(desired_path);
24 if (file.good() && file.is_open()) {
32void Model::init_model() {
33 pinocchio::urdf::buildModel(this->
get_urdf_path(), this->robot_model_);
34 this->robot_data_ = pinocchio::Data(this->robot_model_);
36 std::vector<std::string> frames;
37 for (
auto& f : this->robot_model_.frames) {
38 frames.push_back(f.name);
41 this->frames_ = std::vector<std::string>(frames.begin() + 2, frames.end());
42 this->init_qp_solver();
45bool Model::init_qp_solver() {
47 this->solver_.data()->clearHessianMatrix();
48 this->solver_.data()->clearLinearConstraintsMatrix();
49 this->solver_.clearSolver();
53 this->hessian_ = Eigen::SparseMatrix<double>(nb_joints + 1, nb_joints + 1);
54 this->gradient_ = Eigen::VectorXd::Zero(nb_joints + 1);
55 this->constraint_matrix_ = Eigen::SparseMatrix<double>(3 * nb_joints + 1 + 2, nb_joints + 1);
56 this->lower_bound_constraints_ = Eigen::VectorXd::Zero(3 * nb_joints + 1 + 2);
57 this->upper_bound_constraints_ = Eigen::VectorXd::Zero(3 * nb_joints + 1 + 2);
60 this->hessian_.reserve(nb_joints * nb_joints + 1);
61 this->constraint_matrix_.reserve(5 * nb_joints + 2 * (nb_joints * nb_joints + nb_joints) + 4 * nb_joints + 3);
63 Eigen::VectorXd lower_position_limit = this->robot_model_.lowerPositionLimit;
64 Eigen::VectorXd upper_position_limit = this->robot_model_.upperPositionLimit;
65 Eigen::VectorXd velocity_limit = this->robot_model_.velocityLimit;
68 this->solver_.settings()->setVerbosity(
false);
69 this->solver_.settings()->setWarmStart(
true);
72 for (
unsigned int n = 0; n < nb_joints; ++n) {
74 this->constraint_matrix_.coeffRef(n, n) = 1.0;
76 this->constraint_matrix_.coeffRef(n + nb_joints, n) = 1.0;
77 this->constraint_matrix_.coeffRef(n + nb_joints, nb_joints) = velocity_limit(n);
78 this->upper_bound_constraints_(n + nb_joints) = std::numeric_limits<double>::infinity();
79 this->constraint_matrix_.coeffRef(n + 2 * nb_joints, n) = 1.0;
80 this->constraint_matrix_.coeffRef(n + 2 * nb_joints, nb_joints) = -velocity_limit(n);
81 this->lower_bound_constraints_(n + 2 * nb_joints) = -std::numeric_limits<double>::infinity();
85 this->constraint_matrix_.coeffRef(3 * nb_joints, nb_joints) = 1.0;
86 this->upper_bound_constraints_(3 * nb_joints) = std::numeric_limits<double>::infinity();
88 this->upper_bound_constraints_(3 * nb_joints + 1) = std::numeric_limits<double>::infinity();
89 this->upper_bound_constraints_(3 * nb_joints + 2) = std::numeric_limits<double>::infinity();
92 this->solver_.data()->setNumberOfVariables(
static_cast<int>(nb_joints) + 1);
93 this->solver_.data()->setNumberOfConstraints(this->lower_bound_constraints_.size());
94 if (!this->solver_.data()->setHessianMatrix(this->hessian_)) {
return false; }
95 if (!this->solver_.data()->setGradient(this->gradient_)) {
return false; }
96 if (!this->solver_.data()->setLinearConstraintsMatrix(this->constraint_matrix_)) {
return false; }
97 if (!this->solver_.data()->setLowerBound(this->lower_bound_constraints_)) {
return false; }
98 if (!this->solver_.data()->setUpperBound(this->upper_bound_constraints_)) {
return false; }
100 return this->solver_.initSolver();
103std::vector<unsigned int> Model::get_frame_ids(
const std::vector<std::string>& frames) {
104 std::vector<unsigned int> frame_ids;
105 frame_ids.reserve(frames.size());
107 for (
auto& frame : frames) {
110 frame_ids.push_back(this->robot_model_.frames.size() - 1);
113 if (!this->robot_model_.existFrame(frame)) {
114 throw (exceptions::FrameNotFoundException(frame));
116 frame_ids.push_back(this->robot_model_.getFrameId(frame));
122unsigned int Model::get_frame_id(
const std::string& frame) {
123 return get_frame_ids(std::vector<std::string>{frame}).back();
127 unsigned int frame_id) {
128 if (joint_positions.get_size() != this->get_number_of_joints()) {
129 throw (exceptions::InvalidJointStateSizeException(joint_positions.get_size(), this->get_number_of_joints()));
134 pinocchio::computeFrameJacobian(this->robot_model_,
136 joint_positions.data(),
138 pinocchio::LOCAL_WORLD_ALIGNED,
143 this->robot_model_.frames[frame_id].name,
145 this->get_base_frame());
149 const std::string& frame) {
150 auto frame_id = get_frame_id(frame);
151 return this->compute_jacobian(joint_positions, frame_id);
156 unsigned int frame_id) {
157 if (joint_positions.get_size() != this->get_number_of_joints()) {
160 if (joint_velocities.get_size() != this->get_number_of_joints()) {
161 throw (exceptions::InvalidJointStateSizeException(joint_velocities.get_size(), this->get_number_of_joints()));
165 pinocchio::computeJointJacobiansTimeVariation(this->robot_model_,
167 joint_positions.data(),
168 joint_velocities.data());
169 pinocchio::getFrameJacobianTimeVariation(this->robot_model_,
172 pinocchio::LOCAL_WORLD_ALIGNED,
180 const std::string& frame) {
181 auto frame_id = get_frame_id(frame);
182 return this->compute_jacobian_time_derivative(joint_positions, joint_velocities, frame_id);
187 pinocchio::crba(this->robot_model_, this->robot_data_, joint_positions.data());
189 this->robot_data_.M.triangularView<Eigen::StrictlyLower>() =
190 this->robot_data_.M.transpose().triangularView<Eigen::StrictlyLower>();
191 return this->robot_data_.M;
197 joint_state.get_names(),
198 inertia * joint_state.get_accelerations());
202 return pinocchio::computeCoriolisMatrix(this->robot_model_,
204 joint_state.get_positions(),
205 joint_state.get_velocities());
212 joint_state.get_names(),
213 coriolis_matrix * joint_state.get_velocities());
218 Eigen::VectorXd gravity_torque =
219 pinocchio::computeGeneralizedGravity(this->robot_model_, this->robot_data_, joint_positions.data());
224 unsigned int frame_id) {
225 return this->forward_kinematics(joint_positions, std::vector<unsigned int>{frame_id}).front();
229 const std::vector<unsigned int>& frame_ids) {
230 if (joint_positions.get_size() != this->get_number_of_joints()) {
231 throw (exceptions::InvalidJointStateSizeException(joint_positions.get_size(), this->get_number_of_joints()));
233 std::vector<state_representation::CartesianPose> pose_vector;
234 pinocchio::forwardKinematics(this->robot_model_, this->robot_data_, joint_positions.data());
235 for (
unsigned int id : frame_ids) {
236 if (
id >=
static_cast<unsigned int>(this->robot_model_.nframes)) {
237 throw (exceptions::FrameNotFoundException(std::to_string(
id)));
239 pinocchio::updateFramePlacement(this->robot_model_, this->robot_data_,
id);
240 pinocchio::SE3 pose = this->robot_data_.oMf[id];
241 Eigen::Vector3d translation = pose.translation();
242 Eigen::Quaterniond quaternion;
243 pinocchio::quaternion::assignQuaternion(quaternion, pose.rotation());
247 this->get_base_frame());
248 pose_vector.push_back(frame_pose);
254 const std::string& frame) {
255 std::string actual_frame = frame.empty() ? this->robot_model_.frames.back().name : frame;
256 return this->forward_kinematics(joint_positions, std::vector<std::string>{actual_frame}).front();
260 const std::vector<std::string>& frames) {
261 auto frame_ids = get_frame_ids(frames);
262 return this->forward_kinematics(joint_positions, frame_ids);
266 const double margin) {
267 Eigen::MatrixXd W_b = Eigen::MatrixXd::Identity(this->robot_model_.nq, this->robot_model_.nq);
268 for (
int n = 0; n < this->robot_model_.nq; ++n) {
271 if (joint_positions.data()[n] < this->robot_model_.lowerPositionLimit[n] + margin) {
272 if (joint_positions.data()[n] < this->robot_model_.lowerPositionLimit[n]) {
275 d = (this->robot_model_.lowerPositionLimit[n] + margin - joint_positions.data()[n]) / margin;
276 W_b(n, n) = -2 * d * d * d + 3 * d * d;
278 }
else if (this->robot_model_.upperPositionLimit[n] - margin < joint_positions.data()[n]) {
279 if (this->robot_model_.upperPositionLimit[n] < joint_positions.data()[n]) {
282 d = (joint_positions.data()[n] - (this->robot_model_.upperPositionLimit[n] - margin)) / margin;
283 W_b(n, n) = -2 * d * d * d + 3 * d * d;
292 Eigen::VectorXd Psi(this->robot_model_.nq);
293 Eigen::VectorXd q = joint_positions.data();
294 for (
int i = 0; i < this->robot_model_.nq; ++i) {
296 if (q[i] < this->robot_model_.lowerPositionLimit[i] + margin) {
297 Psi[i] = this->robot_model_.upperPositionLimit[i] - margin
298 - std::max(q[i], this->robot_model_.lowerPositionLimit[i]);
299 }
else if (this->robot_model_.upperPositionLimit[i] - margin < q[i]) {
300 Psi[i] = this->robot_model_.lowerPositionLimit[i] + margin
301 - std::min(q[i], this->robot_model_.upperPositionLimit[i]);
311 const std::string& frame) {
312 std::string actual_frame = frame.empty() ? this->robot_model_.frames.back().name : frame;
313 if (!this->robot_model_.existFrame(actual_frame)) {
317 const std::chrono::nanoseconds dt(
static_cast<int>(1e9));
326 Eigen::MatrixXd JJt(6, 6);
330 Eigen::VectorXd err(6);
331 for (
unsigned int i = 0; i < parameters.max_number_of_iterations; ++i) {
332 err = ((cartesian_pose - this->forward_kinematics(q, actual_frame)) / dt).data();
334 if (err.cwiseAbs().maxCoeff() < parameters.tolerance) {
337 J = this->compute_jacobian(q, actual_frame);
338 W_b = this->cwln_weighted_matrix(q, parameters.margin);
340 psi = parameters.gamma * this->cwln_repulsive_potential_field(q, parameters.margin);
343 JJt.diagonal().array() += parameters.damp;
344 dq.
set_velocities(W_c * psi + parameters.alpha * W_b * (J_b.transpose() * JJt.ldlt().solve(err - J * W_c * psi)));
346 q = this->clamp_in_range(q);
349 err.array().abs().maxCoeff()));
355 const std::string& frame) {
356 Eigen::VectorXd q(pinocchio::neutral(this->robot_model_));
361std::vector<state_representation::CartesianTwist>
363 const std::vector<std::string>& frames) {
364 std::vector<state_representation::CartesianTwist> cartesian_twists(frames.size());
365 for (std::size_t i = 0; i < frames.size(); ++i) {
366 cartesian_twists.at(i) = this->compute_jacobian(joint_state, frames.at(i))
369 return cartesian_twists;
373 const std::string& frame) {
374 return this->
forward_velocity(joint_state, std::vector<std::string>{frame}).front();
377void Model::check_inverse_velocity_arguments(
const std::vector<state_representation::CartesianTwist>& cartesian_twists,
379 const std::vector<std::string>& frames) {
380 if (cartesian_twists.size() != frames.size()) {
381 throw (std::invalid_argument(
"The number of provided twists and frames does not match"));
383 if (joint_positions.get_size() != this->get_number_of_joints()) {
384 throw (exceptions::InvalidJointStateSizeException(joint_positions.get_size(), this->get_number_of_joints()));
386 for (
auto& frame : frames) {
387 if (!this->robot_model_.existFrame(frame)) {
388 throw (exceptions::FrameNotFoundException(frame));
396 const std::vector<std::string>& frames,
397 const double dls_lambda) {
399 this->check_inverse_velocity_arguments(cartesian_twists, joint_positions, frames);
403 Eigen::VectorXd dX(3 * cartesian_twists.size() + 3);
404 Eigen::MatrixXd jacobian(3 * cartesian_twists.size() + 3, nb_joints);
405 for (
unsigned int i = 0; i < cartesian_twists.size() - 1; ++i) {
407 dX.segment<3>(3 * i) = cartesian_twists[i].get_linear_velocity();
408 jacobian.block(3 * i, 0, 3 * i + 3, nb_joints) =
409 this->compute_jacobian(joint_positions, frames.at(i)).
data().block(0, 0, 3, nb_joints);
412 dX.tail(6) = cartesian_twists.back().data();
413 jacobian.bottomRows(6) = this->compute_jacobian(joint_positions, frames.back()).
data();
415 if (dls_lambda == 0.0){
417 joint_positions.get_names(),
418 jacobian.colPivHouseholderQr().solve(dX));
422 if (jacobian.rows() > jacobian.cols()){
423 Eigen::MatrixXd j_prime = jacobian.transpose() * jacobian +
424 dls_lambda * dls_lambda * Eigen::MatrixXd::Identity(jacobian.cols(), jacobian.cols());
426 joint_positions.get_names(),
427 j_prime.colPivHouseholderQr().solve(jacobian.transpose() * dX));
429 Eigen::MatrixXd j_prime = jacobian * jacobian.transpose() +
430 dls_lambda * dls_lambda * Eigen::MatrixXd::Identity(jacobian.rows(), jacobian.rows());
432 joint_positions.get_names(),
433 jacobian.transpose() * j_prime.colPivHouseholderQr().solve(dX));
439 const std::string& frame,
440 const double dls_lambda) {
441 std::string actual_frame = frame.empty() ? this->robot_model_.frames.back().name : frame;
442 return this->
inverse_velocity(std::vector<state_representation::CartesianTwist>({cartesian_twist}),
444 std::vector<std::string>({actual_frame}),
452 const std::vector<std::string>& frames) {
454 using namespace std::chrono;
456 this->check_inverse_velocity_arguments(cartesian_twists, joint_positions, frames);
460 Eigen::VectorXd delta_r(3 * cartesian_twists.size() + 3);
461 Eigen::MatrixXd jacobian(3 * cartesian_twists.size() + 3, nb_joints);
462 for (
unsigned int i = 0; i < cartesian_twists.size() - 1; ++i) {
466 delta_r.segment<3>(3 * i) = displacement.
get_position();
467 jacobian.block(3 * i, 0, 3 * i + 3, nb_joints) =
468 this->compute_jacobian(joint_positions, frames.at(i)).
data().block(0, 0, 3, nb_joints);
472 delta_r.segment<3>(3 * (cartesian_twists.size() - 1)) = full_displacement.
get_position();
474 jacobian.bottomRows(6) = this->compute_jacobian(joint_positions, frames.back()).
data();
476 Eigen::MatrixXd hessian_matrix = jacobian.transpose() * jacobian;
478 std::vector<Eigen::Triplet<double>> coefficients;
479 for (
unsigned int i = 0; i < nb_joints; ++i) {
480 for (
unsigned int j = 0; j < nb_joints; ++j) {
481 coefficients.emplace_back(Eigen::Triplet<double>(i, j, hessian_matrix(i, j)));
484 coefficients.emplace_back(Eigen::Triplet<double>(nb_joints, nb_joints, parameters.alpha));
485 this->hessian_.setFromTriplets(coefficients.begin(), coefficients.end());
487 this->gradient_.head(nb_joints) = -parameters.proportional_gain * delta_r.transpose() * jacobian;
489 this->lower_bound_constraints_(3 * nb_joints) = duration_cast<duration<float>>(parameters.dt).count();
491 Eigen::VectorXd lower_position_limit = this->robot_model_.lowerPositionLimit;
492 Eigen::VectorXd upper_position_limit = this->robot_model_.upperPositionLimit;
493 for (
unsigned int n = 0; n < nb_joints; ++n) {
494 this->lower_bound_constraints_(n) = lower_position_limit(n) - joint_positions.data()(n);
495 this->upper_bound_constraints_(n) = upper_position_limit(n) - joint_positions.data()(n);
498 this->constraint_matrix_.coeffRef(3 * nb_joints + 1, nb_joints) = parameters.linear_velocity_limit;
499 this->constraint_matrix_.coeffRef(3 * nb_joints + 2, nb_joints) = parameters.angular_velocity_limit;
500 this->lower_bound_constraints_(3 * nb_joints + 1) = full_displacement.
get_position().norm();
501 this->lower_bound_constraints_(3 * nb_joints + 2) = full_displacement.
get_orientation().vec().norm();
504 this->solver_.updateHessianMatrix(this->hessian_);
505 this->solver_.updateGradient(this->gradient_);
506 this->solver_.updateBounds(this->lower_bound_constraints_, this->upper_bound_constraints_);
507 this->solver_.updateLinearConstraintsMatrix(this->constraint_matrix_);
509 this->solver_.solve();
512 joint_positions.get_names(),
513 this->solver_.getSolution().head(nb_joints));
514 double dt = this->solver_.getSolution().tail(1)(0);
521 const std::string& frame) {
522 std::string actual_frame = frame.empty() ? this->robot_model_.frames.back().name : frame;
523 return this->
inverse_velocity(std::vector<state_representation::CartesianTwist>({cartesian_twist}),
526 std::vector<std::string>({actual_frame}));
530 std::cout <<
"hessian:" << std::endl;
531 std::cout << this->hessian_ << std::endl;
533 std::cout <<
"gradient:" << std::endl;
534 std::cout << this->gradient_ << std::endl;
536 for (
unsigned int i = 0; i < this->constraint_matrix_.rows(); ++i) {
537 std::cout << this->lower_bound_constraints_(i);
538 std::cout <<
" < | ";
539 for (
unsigned int j = 0; j < this->constraint_matrix_.cols(); ++j) {
540 std::cout << this->constraint_matrix_.coeffRef(i, j) <<
" | ";
543 std::cout << this->upper_bound_constraints_(i) << std::endl;
547bool Model::in_range(
const Eigen::VectorXd& vector,
548 const Eigen::VectorXd& lower_limits,
549 const Eigen::VectorXd& upper_limits) {
550 return ((vector.array() >= lower_limits.array()).all() && (vector.array() <= upper_limits.array()).all());
554 return this->in_range(joint_positions.get_positions(),
555 this->robot_model_.lowerPositionLimit,
556 this->robot_model_.upperPositionLimit);
560 return this->in_range(joint_velocities.get_velocities(),
561 -this->robot_model_.velocityLimit,
562 this->robot_model_.velocityLimit);
566 return this->in_range(joint_torques.get_torques(), -this->robot_model_.effortLimit, this->robot_model_.effortLimit);
575Eigen::VectorXd Model::clamp_in_range(
const Eigen::VectorXd& vector,
576 const Eigen::VectorXd& lower_limits,
577 const Eigen::VectorXd& upper_limits) {
578 return lower_limits.cwiseMax(upper_limits.cwiseMin(vector));
583 joint_state_clamped.
set_positions(this->clamp_in_range(joint_state.get_positions(),
584 this->robot_model_.lowerPositionLimit,
585 this->robot_model_.upperPositionLimit));
586 joint_state_clamped.
set_velocities(this->clamp_in_range(joint_state.get_velocities(),
587 -this->robot_model_.velocityLimit,
588 this->robot_model_.velocityLimit));
589 joint_state_clamped.
set_torques(this->clamp_in_range(joint_state.get_torques(),
590 -this->robot_model_.effortLimit,
591 this->robot_model_.effortLimit));
592 return joint_state_clamped;
The Model class is a wrapper around pinocchio dynamic computation library with state_representation e...
unsigned int get_number_of_joints() const
Getter of the number of joints.
std::vector< state_representation::CartesianTwist > forward_velocity(const state_representation::JointState &joint_state, const std::vector< std::string > &frames)
Compute the forward velocity kinematics, i.e. the twist of certain frames from the joint states.
Eigen::MatrixXd compute_coriolis_matrix(const state_representation::JointState &joint_state)
Compute the Coriolis matrix from a given joint state.
state_representation::JointTorques compute_inertia_torques(const state_representation::JointState &joint_state)
Compute the Inertia torques, i.e the inertia matrix multiplied by the joint accelerations....
state_representation::JointVelocities inverse_velocity(const std::vector< state_representation::CartesianTwist > &cartesian_twists, const state_representation::JointPositions &joint_positions, const std::vector< std::string > &frames, const double dls_lambda=0.0)
Compute the inverse velocity kinematics, i.e. joint velocities from the velocities of the frames in p...
const std::string & get_base_frame() const
Getter of the base frame of the robot.
state_representation::JointTorques compute_gravity_torques(const state_representation::JointPositions &joint_positions)
Compute the gravity torques.
Eigen::MatrixXd compute_inertia_matrix(const state_representation::JointPositions &joint_positions)
Compute the Inertia matrix from a given joint positions.
Model(const std::string &robot_name, const std::string &urdf_path)
Constructor with robot name and path to URDF file.
const std::string & get_urdf_path() const
Getter of the URDF path.
state_representation::JointPositions inverse_kinematics(const state_representation::CartesianPose &cartesian_pose, const InverseKinematicsParameters ¶meters=InverseKinematicsParameters(), const std::string &frame="")
Compute the inverse kinematics, i.e. joint positions from the pose of the end-effector in an iterativ...
state_representation::JointTorques compute_coriolis_torques(const state_representation::JointState &joint_state)
Compute the Coriolis torques, i.e. the Coriolis matrix multiplied by the joint velocities and express...
void print_qp_problem()
Helper function to print the qp_problem (for debugging)
const std::string & get_robot_name() const
Getter of the robot name.
std::vector< std::string > get_joint_frames() const
Getter of the joint frames from the model.
static bool create_urdf_from_string(const std::string &urdf_string, const std::string &desired_path)
Creates a URDF file with desired path and name from a string (possibly the robot description string f...
Class to define Cartesian pose in Cartesian space as 3D position and quaternion based orientation.
const Eigen::Vector3d & get_position() const
Getter of the position attribute.
const Eigen::Quaterniond & get_orientation() const
Getter of the orientation attribute.
Class to define twist in Cartesian space as 3D linear and angular velocity vectors.
Class to define a robot Jacobian matrix.
const Eigen::MatrixXd & data() const
Getter of the data attribute.
Eigen::MatrixXd transpose() const
Return the transpose of the Jacobian matrix.
Class to define positions of the joints.
Class to define a state in joint space.
void set_positions(const Eigen::VectorXd &positions)
Setter of the positions attribute.
void set_torques(const Eigen::VectorXd &torques)
Setter of the torques attribute.
const std::vector< std::string > & get_names() const
Getter of the names attribute.
void set_velocities(const Eigen::VectorXd &velocities)
Setter of the velocities attribute.
Class to define torques of the joints.
Class to define velocities of the joints.
Robot kinematics and dynamics.
Core state variables and objects.
parameters for the inverse kinematics function
parameters for the inverse velocity kinematics function