28
28
#include " pinocchio/algorithm/frames.hpp"
29
29
#include " pinocchio/algorithm/jacobian.hpp"
30
30
#include " pinocchio/algorithm/kinematics.hpp"
31
+ #include " pseudoinverse.hpp"
31
32
#include " tf2_eigen/tf2_eigen.hpp"
32
33
33
34
namespace ik_solvers
@@ -173,7 +174,7 @@ namespace
173
174
auto compute_nullspace (const Eigen::MatrixXd & augmented_jacobian) -> Eigen::MatrixXd
174
175
{
175
176
const Eigen::MatrixXd eye = Eigen::MatrixXd::Identity (augmented_jacobian.cols (), augmented_jacobian.cols ());
176
- return eye - (augmented_jacobian. completeOrthogonalDecomposition (). pseudoInverse ( ) * augmented_jacobian);
177
+ return eye - (pinv::least_squares (augmented_jacobian ) * augmented_jacobian);
177
178
}
178
179
179
180
// / Construct the augmented Jacobian matrix from a list of Jacobian matrices.
@@ -198,42 +199,37 @@ auto construct_augmented_jacobian(const std::vector<Eigen::MatrixXd> & jacobians
198
199
return augmented_jacobian;
199
200
}
200
201
201
- // / Closed-loop task priority IK using the damped pseudoinverse .
202
+ // / Closed-loop task priority IK.
202
203
auto tpik (const hierarchy::ConstraintSet & tasks, size_t nv, double damping)
203
204
-> std::expected<Eigen::VectorXd, SolverError>
204
205
{
205
206
if (tasks.empty ()) {
206
207
return std::unexpected (SolverError::NO_SOLUTION);
207
208
}
208
209
209
- Eigen::VectorXd vel = Eigen::VectorXd::Zero (nv);
210
210
std::vector<Eigen::MatrixXd> jacobians;
211
211
jacobians.reserve (tasks.size ());
212
+
213
+ Eigen::VectorXd vel = Eigen::VectorXd::Zero (nv);
212
214
Eigen::MatrixXd nullspace = Eigen::MatrixXd::Identity (nv, nv);
213
215
214
216
for (const auto & task : tasks) {
215
- // damped pseudoinverse of the Jacobian
216
- const Eigen::MatrixXd J = task->jacobian ();
217
- const auto eye = Eigen::MatrixXd::Identity (J.rows (), J.rows ());
218
- const Eigen::MatrixXd J_inv = J.transpose () * (J * J.transpose () + damping * eye).inverse ();
219
-
220
- // closed-loop task priority IK
217
+ const Eigen::MatrixXd J_inv = pinv::damped_least_squares (task->jacobian (), damping);
221
218
vel += nullspace * J_inv * (task->gain () * task->error ());
222
219
223
- jacobians.push_back (J );
220
+ jacobians.push_back (task-> jacobian () );
224
221
nullspace = compute_nullspace (construct_augmented_jacobian (jacobians));
225
222
}
226
223
227
224
return vel;
228
225
}
229
226
230
- // / Check if the solution is feasible with respect to the set constraints .
227
+ // / Check if the solution is feasible.
231
228
auto is_feasible (const hierarchy::ConstraintSet & constraints, const Eigen::VectorXd & solution) -> bool
232
229
{
233
230
for (const auto & constraint : constraints) {
234
231
auto set_task = std::dynamic_pointer_cast<hierarchy::SetConstraint>(constraint);
235
232
const double pred = (constraint->jacobian () * solution).value ();
236
-
237
233
if (!((set_task->primal () > set_task->lower_threshold () && pred < 0 ) ||
238
234
(set_task->primal () < set_task->upper_threshold () && pred > 0 ))) {
239
235
return false ;
@@ -253,39 +249,32 @@ auto search_solutions(
253
249
return std::unexpected (SolverError::NO_SOLUTION);
254
250
}
255
251
256
- Eigen::VectorXd solution;
257
-
252
+ // there is only one solution if there are no set tasks
258
253
if (set_tasks.empty ()) {
259
- const auto out = tpik (hierarchies.front (), nv, damping);
260
- if (!out.has_value ()) {
261
- return std::unexpected (out.error ());
262
- }
263
- solution = out.value ();
264
- } else {
265
- std::vector<Eigen::VectorXd> solutions;
266
- solutions.reserve (hierarchies.size ());
267
-
268
- for (const auto & tasks : hierarchies) {
269
- const auto out = tpik (tasks, nv, damping);
270
- if (!out.has_value ()) {
271
- continue ;
272
- }
254
+ return tpik (hierarchies.front (), nv, damping);
255
+ }
273
256
274
- const Eigen::VectorXd & current_solution = out.value ();
275
- if (is_feasible (set_tasks, current_solution)) {
276
- solutions.push_back (current_solution);
277
- }
257
+ std::vector<Eigen::VectorXd> solutions;
258
+ solutions.reserve (hierarchies.size ());
259
+
260
+ for (const auto & tasks : hierarchies) {
261
+ const auto out = tpik (tasks, nv, damping);
262
+ if (!out.has_value ()) {
263
+ continue ;
278
264
}
279
265
280
- if (solutions.empty ()) {
281
- return std::unexpected (SolverError::NO_SOLUTION);
266
+ const Eigen::VectorXd & current_solution = out.value ();
267
+ if (is_feasible (set_tasks, current_solution)) {
268
+ solutions.push_back (current_solution);
282
269
}
270
+ }
283
271
284
- // Choose the solution with the smallest norm
285
- solution = * std::ranges::min_element (solutions, {}, []( const auto & a) { return a. norm (); } );
272
+ if (solutions. empty ()) {
273
+ return std::unexpected (SolverError::NO_SOLUTION );
286
274
}
287
275
288
- return solution;
276
+ // Choose the solution with the smallest norm
277
+ return *std::ranges::min_element (solutions, {}, [](const auto & a) { return a.norm (); });
289
278
}
290
279
291
280
auto pinocchio_to_eigen (const pinocchio::SE3 & pose) -> Eigen::Affine3d
@@ -354,7 +343,7 @@ auto TaskPriorityIKSolver::solve_ik(const Eigen::Affine3d & target_pose, const E
354
343
const double activation = params_.joint_limit_task .constrained_joints_map [joint_name].activation_threshold ;
355
344
const double joint_limit_gain = params_.joint_limit_task .constrained_joints_map [joint_name].gain ;
356
345
357
- // insert the joint constraint
346
+ // insert the joint limit constraint
358
347
hierarchy_.insert (std::make_shared<hierarchy::JointConstraint>(
359
348
model_, primal, ub, lb, tol, activation, joint_name, joint_limit_gain));
360
349
}
0 commit comments