GCC Code Coverage Report


Directory: ./
File: include/crocoddyl/multibody/actuations/floating-base-thrusters.hpp
Date: 2025-02-24 23:41:29
Exec Total Coverage
Lines: 81 120 67.5%
Branches: 68 209 32.5%

Line Branch Exec Source
1 ///////////////////////////////////////////////////////////////////////////////
2 // BSD 3-Clause License
3 //
4 // Copyright (C) 2014-2024, Heriot-Watt University
5 // Copyright note valid unless otherwise stated in individual files.
6 // All rights reserved.
7 ///////////////////////////////////////////////////////////////////////////////
8
9 #ifndef CROCODDYL_MULTIBODY_ACTUATIONS_PROPELLERS_HPP_
10 #define CROCODDYL_MULTIBODY_ACTUATIONS_PROPELLERS_HPP_
11
12 #include <pinocchio/multibody/fwd.hpp>
13 #include <pinocchio/spatial/se3.hpp>
14 #include <vector>
15
16 #include "crocoddyl/core/actuation-base.hpp"
17 #include "crocoddyl/core/utils/exception.hpp"
18 #include "crocoddyl/multibody/states/multibody.hpp"
19
20 namespace crocoddyl {
21
22 enum ThrusterType { CW = 0, CCW };
23
24 template <typename _Scalar>
25 struct ThrusterTpl {
26 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
27
28 typedef _Scalar Scalar;
29 typedef pinocchio::SE3Tpl<Scalar> SE3;
30 typedef ThrusterTpl<Scalar> Thruster;
31
32 /**
33 * @brief Initialize the thruster in a give pose from the root joint.
34 *
35 * @param pose[in] Pose from root joint
36 * @param ctorque[in] Coefficient of generated torque per thrust
37 * @param type[in] Type of thruster (clockwise or counterclockwise,
38 * default clockwise)
39 * @param[in] min_thrust[in] Minimum thrust (default 0.)
40 * @param[in] max_thrust[in] Maximum thrust (default inf number))
41 */
42 6840 ThrusterTpl(const SE3& pose, const Scalar ctorque,
43 const ThrusterType type = CW,
44 const Scalar min_thrust = Scalar(0.),
45 const Scalar max_thrust = std::numeric_limits<Scalar>::infinity())
46 6840 : pose(pose),
47 6840 ctorque(ctorque),
48 6840 type(type),
49 6840 min_thrust(min_thrust),
50 6840 max_thrust(max_thrust) {}
51
52 /**
53 * @brief Initialize the thruster in a pose in the origin of the root joint.
54 *
55 * @param pose[in] Pose from root joint
56 * @param ctorque[in] Coefficient of generated torque per thrust
57 * @param type[in] Type of thruster (clockwise or counterclockwise,
58 * default clockwise)
59 * @param[in] min_thrust[in] Minimum thrust (default 0.)
60 * @param[in] max_thrust[in] Maximum thrust (default inf number))
61 */
62 ThrusterTpl(const Scalar ctorque, const ThrusterType type = CW,
63 const Scalar min_thrust = Scalar(0.),
64 const Scalar max_thrust = std::numeric_limits<Scalar>::infinity())
65 : pose(SE3::Identity()),
66 ctorque(ctorque),
67 type(type),
68 min_thrust(min_thrust),
69 max_thrust(max_thrust) {}
70 12793 ThrusterTpl(const ThrusterTpl<Scalar>& clone)
71 12793 : pose(clone.pose),
72 12793 ctorque(clone.ctorque),
73 12793 type(clone.type),
74 12793 min_thrust(clone.min_thrust),
75 12793 max_thrust(clone.max_thrust) {}
76
77 ThrusterTpl& operator=(const ThrusterTpl<Scalar>& other) {
78 if (this != &other) {
79 pose = other.pose;
80 ctorque = other.ctorque;
81 type = other.type;
82 min_thrust = other.min_thrust;
83 max_thrust = other.max_thrust;
84 }
85 return *this;
86 }
87
88 template <typename OtherScalar>
89 bool operator==(const ThrusterTpl<OtherScalar>& other) const {
90 return (pose == other.pose && ctorque == other.ctorque &&
91 type == other.type && min_thrust == other.min_thrust &&
92 max_thrust == other.max_thrust);
93 }
94
95 friend std::ostream& operator<<(std::ostream& os,
96 const ThrusterTpl<Scalar>& X) {
97 os << " pose:" << std::endl
98 << X.pose << " ctorque: " << X.ctorque << std::endl
99 << " type: " << X.type << std::endl
100 << "min_thrust: " << X.min_thrust << std::endl
101 << "max_thrust: " << X.max_thrust << std::endl;
102 return os;
103 }
104
105 SE3 pose; //!< Thruster pose
106 Scalar ctorque; //!< Coefficient of generated torque per thrust
107 ThrusterType type; //!< Type of thruster (CW and CCW for clockwise and
108 //!< counterclockwise, respectively)
109 Scalar min_thrust; //!< Minimum thrust
110 Scalar max_thrust; //!< Minimum thrust
111 };
112
113 /**
114 * @brief Actuation models for floating base systems actuated with thrusters
115 *
116 * This actuation model models floating base robots equipped with thrusters,
117 * e.g., multicopters or marine robots equipped with manipulators. It control
118 * inputs are the thrusters' thrust (i.e., forces) and joint torques.
119 *
120 * Both actuation and Jacobians are computed analytically by `calc` and
121 * `calcDiff`, respectively.
122 *
123 * We assume the robot velocity to zero for easily related square thruster
124 * velocities with thrust and torque generated. This approach is similarly
125 * implemented in M. Geisert and N. Mansard, "Trajectory generation for
126 * quadrotor based systems using numerical optimal control", (ICRA). See Section
127 * III.C.
128 *
129 * \sa `ActuationModelAbstractTpl`, `calc()`, `calcDiff()`, `createData()`
130 */
131 template <typename _Scalar>
132 class ActuationModelFloatingBaseThrustersTpl
133 : public ActuationModelAbstractTpl<_Scalar> {
134 public:
135 typedef _Scalar Scalar;
136 typedef MathBaseTpl<Scalar> MathBase;
137 typedef ActuationModelAbstractTpl<Scalar> Base;
138 typedef ActuationDataAbstractTpl<Scalar> Data;
139 typedef StateMultibodyTpl<Scalar> StateMultibody;
140 typedef ThrusterTpl<Scalar> Thruster;
141 typedef typename MathBase::Vector3s Vector3s;
142 typedef typename MathBase::VectorXs VectorXs;
143 typedef typename MathBase::MatrixXs MatrixXs;
144
145 /**
146 * @brief Initialize the floating base actuation model equipped with
147 * thrusters
148 *
149 * @param[in] state State of the dynamical system
150 * @param[in] thrusters Vector of thrusters
151 */
152 198 ActuationModelFloatingBaseThrustersTpl(std::shared_ptr<StateMultibody> state,
153 const std::vector<Thruster>& thrusters)
154 : Base(state,
155 198 state->get_nv() -
156 198 state->get_pinocchio()
157 198 ->joints[(
158
3/6
✓ Branch 5 taken 198 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 198 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 198 times.
✗ Branch 11 not taken.
396 state->get_pinocchio()->existJointName("root_joint")
159
3/10
✓ Branch 5 taken 198 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 198 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 198 times.
✗ Branch 12 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
396 ? state->get_pinocchio()->getJointId("root_joint")
160 : 0)]
161
1/2
✓ Branch 1 taken 198 times.
✗ Branch 2 not taken.
198 .nv() +
162 198 thrusters.size()),
163 198 thrusters_(thrusters),
164 198 n_thrusters_(thrusters.size()),
165
1/2
✓ Branch 3 taken 198 times.
✗ Branch 4 not taken.
198 W_thrust_(state_->get_nv(), nu_),
166
4/8
✓ Branch 4 taken 198 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 198 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 198 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 198 times.
✗ Branch 14 not taken.
792 update_data_(true) {
167
3/6
✓ Branch 5 taken 198 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 198 times.
✗ Branch 9 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 198 times.
198 if (!state->get_pinocchio()->existJointName("root_joint")) {
168 throw_pretty(
169 "Invalid argument: "
170 << "the first joint has to be a root one (e.g., free-flyer joint)");
171 }
172 // Update the joint actuation part
173
1/2
✓ Branch 1 taken 198 times.
✗ Branch 2 not taken.
198 W_thrust_.setZero();
174
2/2
✓ Branch 0 taken 92 times.
✓ Branch 1 taken 106 times.
198 if (nu_ > n_thrusters_) {
175
1/2
✓ Branch 1 taken 92 times.
✗ Branch 2 not taken.
92 W_thrust_.bottomRightCorner(nu_ - n_thrusters_, nu_ - n_thrusters_)
176 .diagonal()
177
2/4
✓ Branch 1 taken 92 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 92 times.
✗ Branch 5 not taken.
92 .setOnes();
178 }
179 // Update the floating base actuation part
180
1/2
✓ Branch 1 taken 198 times.
✗ Branch 2 not taken.
198 set_thrusters(thrusters_);
181 198 }
182 400 virtual ~ActuationModelFloatingBaseThrustersTpl() {}
183
184 /**
185 * @brief Compute the actuation signal and actuation set from its thrust
186 * and joint torque inputs \f$\mathbf{u}\in\mathbb{R}^{nu}\f$
187 *
188 * @param[in] data Floating base thrusters actuation data
189 * @param[in] x State point \f$\mathbf{x}\in\mathbb{R}^{ndx}\f$
190 * @param[in] u Joint-torque input \f$\mathbf{u}\in\mathbb{R}^{nu}\f$
191 */
192 9310 virtual void calc(const std::shared_ptr<Data>& data,
193 const Eigen::Ref<const VectorXs>&,
194 const Eigen::Ref<const VectorXs>& u) {
195
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9310 times.
9310 if (static_cast<std::size_t>(u.size()) != nu_) {
196 throw_pretty(
197 "Invalid argument: " << "u has wrong dimension (it should be " +
198 std::to_string(nu_) + ")");
199 }
200
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9310 times.
9310 if (update_data_) {
201 updateData(data);
202 }
203
2/4
✓ Branch 4 taken 9310 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9310 times.
✗ Branch 8 not taken.
9310 data->tau.noalias() = data->dtau_du * u;
204 9310 }
205
206 /**
207 * @brief Compute the Jacobians of the floating base thruster actuation
208 * function
209 *
210 * @param[in] data Floating base thrusters actuation data
211 * @param[in] x State point \f$\mathbf{x}\in\mathbb{R}^{ndx}\f$
212 * @param[in] u Joint-torque input \f$\mathbf{u}\in\mathbb{R}^{nu}\f$
213 */
214 #ifndef NDEBUG
215 1606 virtual void calcDiff(const std::shared_ptr<Data>& data,
216 const Eigen::Ref<const VectorXs>&,
217 const Eigen::Ref<const VectorXs>&) {
218 #else
219 virtual void calcDiff(const std::shared_ptr<Data>&,
220 const Eigen::Ref<const VectorXs>&,
221 const Eigen::Ref<const VectorXs>&) {
222 #endif
223 // The derivatives has constant values which were set in createData.
224
2/12
✓ Branch 4 taken 1606 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 1606 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
1606 assert_pretty(MatrixXs(data->dtau_du).isApprox(W_thrust_),
225 "dtau_du has wrong value");
226 1606 }
227
228 2313 virtual void commands(const std::shared_ptr<Data>& data,
229 const Eigen::Ref<const VectorXs>&,
230 const Eigen::Ref<const VectorXs>& tau) {
231
2/4
✓ Branch 4 taken 2313 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2313 times.
✗ Branch 8 not taken.
2313 data->u.noalias() = data->Mtau * tau;
232 2313 }
233
234 #ifndef NDEBUG
235 629 virtual void torqueTransform(const std::shared_ptr<Data>& data,
236 const Eigen::Ref<const VectorXs>&,
237 const Eigen::Ref<const VectorXs>&) {
238 #else
239 virtual void torqueTransform(const std::shared_ptr<Data>&,
240 const Eigen::Ref<const VectorXs>&,
241 const Eigen::Ref<const VectorXs>&) {
242 #endif
243 // The torque transform has constant values which were set in createData.
244
2/12
✓ Branch 4 taken 629 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 629 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
629 assert_pretty(MatrixXs(data->Mtau).isApprox(Mtau_), "Mtau has wrong value");
245 629 }
246
247 /**
248 * @brief Create the floating base thruster actuation data
249 *
250 * @return the actuation data
251 */
252 6522 virtual std::shared_ptr<Data> createData() {
253
1/2
✓ Branch 1 taken 6522 times.
✗ Branch 2 not taken.
6522 std::shared_ptr<Data> data =
254 6522 std::allocate_shared<Data>(Eigen::aligned_allocator<Data>(), this);
255
1/2
✓ Branch 1 taken 6522 times.
✗ Branch 2 not taken.
6522 updateData(data);
256 6522 return data;
257 }
258
259 /**
260 * @brief Return the vector of thrusters
261 */
262 const std::vector<Thruster>& get_thrusters() const { return thrusters_; }
263
264 /**
265 * @brief Return the number of thrusters
266 */
267 std::size_t get_nthrusters() const { return n_thrusters_; }
268
269 /**
270 * @brief Modify the vector of thrusters
271 *
272 * Since we don't want to allocate data, we request to pass the same
273 * number of thrusters.
274 *
275 * @param[in] thrusters Vector of thrusters
276 */
277 198 void set_thrusters(const std::vector<Thruster>& thrusters) {
278
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 198 times.
198 if (static_cast<std::size_t>(thrusters.size()) != n_thrusters_) {
279 throw_pretty("Invalid argument: "
280 << "the number of thrusters is wrong (it should be " +
281 std::to_string(n_thrusters_) + ")");
282 }
283 198 thrusters_ = thrusters;
284 // Update the mapping matrix from thrusters thrust to body force/moments
285
2/2
✓ Branch 0 taken 792 times.
✓ Branch 1 taken 198 times.
990 for (std::size_t i = 0; i < n_thrusters_; ++i) {
286 792 const Thruster& p = thrusters_[i];
287
4/8
✓ Branch 1 taken 792 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 792 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 792 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 792 times.
✗ Branch 11 not taken.
792 const Vector3s& f_z = p.pose.rotation() * Vector3s::UnitZ();
288
3/6
✓ Branch 1 taken 792 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 792 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 792 times.
✗ Branch 8 not taken.
792 W_thrust_.template topRows<3>().col(i) += f_z;
289
4/8
✓ Branch 1 taken 792 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 792 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 792 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 792 times.
✗ Branch 11 not taken.
792 W_thrust_.template middleRows<3>(3).col(i).noalias() +=
290
2/4
✓ Branch 1 taken 792 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 792 times.
✗ Branch 5 not taken.
792 p.pose.translation().cross(f_z);
291
2/3
✓ Branch 0 taken 396 times.
✓ Branch 1 taken 396 times.
✗ Branch 2 not taken.
792 switch (p.type) {
292 396 case CW:
293
4/8
✓ Branch 1 taken 396 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 396 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 396 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 396 times.
✗ Branch 11 not taken.
396 W_thrust_.template middleRows<3>(3).col(i) += p.ctorque * f_z;
294 396 break;
295 396 case CCW:
296
4/8
✓ Branch 1 taken 396 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 396 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 396 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 396 times.
✗ Branch 11 not taken.
396 W_thrust_.template middleRows<3>(3).col(i) -= p.ctorque * f_z;
297 396 break;
298 }
299 }
300 // Compute the torque transform matrix from generalized torques to joint
301 // torque inputs
302
2/4
✓ Branch 2 taken 198 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 198 times.
✗ Branch 6 not taken.
198 Mtau_ = pseudoInverse(MatrixXs(W_thrust_));
303
2/4
✓ Branch 2 taken 198 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 198 times.
✗ Branch 6 not taken.
198 S_.noalias() = W_thrust_ * Mtau_;
304 198 update_data_ = true;
305 198 }
306
307 const MatrixXs& get_Wthrust() const { return W_thrust_; }
308
309 const MatrixXs& get_S() const { return S_; }
310
311 void print(std::ostream& os) const {
312 os << "ActuationModelFloatingBaseThrusters {nu=" << nu_
313 << ", nthrusters=" << n_thrusters_ << ", thrusters=" << std::endl;
314 for (std::size_t i = 0; i < n_thrusters_; ++i) {
315 os << std::to_string(i) << ": " << thrusters_[i];
316 }
317 os << "}";
318 }
319
320 protected:
321 std::vector<Thruster> thrusters_; //!< Vector of thrusters
322 std::size_t n_thrusters_; //!< Number of thrusters
323 MatrixXs W_thrust_; //!< Matrix from thrusters thrusts to body wrench
324 MatrixXs Mtau_; //!< Constaint torque transform from generalized torques to
325 //!< joint torque inputs
326 MatrixXs S_; //!< Selection matrix for under-actuation part
327
328 bool update_data_;
329 using Base::nu_;
330 using Base::state_;
331
332 private:
333 6522 void updateData(const std::shared_ptr<Data>& data) {
334 6522 data->dtau_du = W_thrust_;
335 6522 data->Mtau = Mtau_;
336 6522 const std::size_t nv = state_->get_nv();
337
2/2
✓ Branch 0 taken 41284 times.
✓ Branch 1 taken 6522 times.
47806 for (std::size_t k = 0; k < nv; ++k) {
338
2/2
✓ Branch 2 taken 13044 times.
✓ Branch 3 taken 28240 times.
41284 if (fabs(S_(k, k)) < std::numeric_limits<Scalar>::epsilon()) {
339 13044 data->tau_set[k] = false;
340 } else {
341 28240 data->tau_set[k] = true;
342 }
343 }
344 6522 update_data_ = false;
345 6522 }
346 };
347
348 } // namespace crocoddyl
349
350 #endif // CROCODDYL_MULTIBODY_ACTUATIONS_PROPELLERS_HPP_
351