GCC Code Coverage Report


Directory: ./
File: include/crocoddyl/core/solver-base.hpp
Date: 2025-03-26 19:23:43
Exec Total Coverage
Lines: 2 2 100.0%
Branches: 0 0 -%

Line Branch Exec Source
1 ///////////////////////////////////////////////////////////////////////////////
2 // BSD 3-Clause License
3 //
4 // Copyright (C) 2019-2025, LAAS-CNRS, University of Edinburgh,
5 // Heriot-Watt University
6 // Copyright note valid unless otherwise stated in individual files.
7 // All rights reserved.
8 ///////////////////////////////////////////////////////////////////////////////
9
10 #ifndef CROCODDYL_CORE_SOLVER_BASE_HPP_
11 #define CROCODDYL_CORE_SOLVER_BASE_HPP_
12
13 #include <vector>
14
15 #include "crocoddyl/core/fwd.hpp"
16 #include "crocoddyl/core/optctrl/shooting.hpp"
17 #include "crocoddyl/core/utils/stop-watch.hpp"
18
19 namespace crocoddyl {
20
21 class CallbackAbstract; // forward declaration
22 static std::vector<Eigen::VectorXd> DEFAULT_VECTOR;
23
24 enum FeasibilityNorm { LInf = 0, L1 };
25
26 /**
27 * @brief Abstract class for optimal control solvers
28 *
29 * A solver resolves an optimal control solver of the form
30 * \f{eqnarray*}{
31 * \begin{Bmatrix}
32 * \mathbf{x}^*_0,\cdots,\mathbf{x}^*_{T} \\
33 * \mathbf{u}^*_0,\cdots,\mathbf{u}^*_{T-1}
34 * \end{Bmatrix} =
35 * \arg\min_{\mathbf{x}_s,\mathbf{u}_s} && l_T (\mathbf{x}_T) + \sum_{k=0}^{T-1}
36 * l_k(\mathbf{x}_t,\mathbf{u}_t) \\
37 * \operatorname{subject}\,\operatorname{to} && \mathbf{x}_0 =
38 * \mathbf{\tilde{x}}_0\\
39 * && \mathbf{x}_{k+1} = \mathbf{f}_k(\mathbf{x}_k,\mathbf{u}_k)\\
40 * && \mathbf{x}_k\in\mathcal{X}, \mathbf{u}_k\in\mathcal{U}
41 * \f}
42 * where \f$l_T(\mathbf{x}_T)\f$, \f$l_k(\mathbf{x}_t,\mathbf{u}_t)\f$ are the
43 * terminal and running cost functions, respectively,
44 * \f$\mathbf{f}_k(\mathbf{x}_k,\mathbf{u}_k)\f$ describes evolution of the
45 * system, and state and control admissible sets are defined by
46 * \f$\mathbf{x}_k\in\mathcal{X}\f$, \f$\mathbf{u}_k\in\mathcal{U}\f$. An action
47 * model, defined in the shooting problem, describes each node \f$k\f$. Inside
48 * the action model, we specialize the cost functions, the system evolution and
49 * the admissible sets.
50 *
51 * The main routines are `computeDirection()` and `tryStep()`. The former finds
52 * a search direction and typically computes the derivatives of each action
53 * model. The latter rollout the dynamics and cost (i.e., the action) to try the
54 * search direction found by `computeDirection`. Both functions used the current
55 * guess defined by `setCandidate()`. Finally, `solve()` function is used to
56 * define when the search direction and length are computed in each iterate. It
57 * also describes the globalization strategy (i.e., regularization) of the
58 * numerical optimization.
59 *
60 * \sa `solve()`, `computeDirection()`, `tryStep()`, `stoppingCriteria()`
61 */
62 class SolverAbstract {
63 public:
64 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
65
66 /**
67 * @brief Initialize the solver
68 *
69 * @param[in] problem shooting problem
70 */
71 explicit SolverAbstract(std::shared_ptr<ShootingProblem> problem);
72 virtual ~SolverAbstract();
73
74 /**
75 * @brief Compute the optimal trajectory \f$\mathbf{x}^*_s,\mathbf{u}^*_s\f$
76 * as lists of \f$T+1\f$ and \f$T\f$ terms
77 *
78 * From an initial guess \p init_xs, \p init_us (feasible or not), iterate
79 * over `computeDirection()` and `tryStep()` until `stoppingCriteria()` is
80 * below threshold. It also describes the globalization strategy used during
81 * the numerical optimization.
82 *
83 * @param[in] init_xs initial guess for state trajectory with \f$T+1\f$
84 * elements (default [])
85 * @param[in] init_us initial guess for control trajectory with \f$T\f$
86 * elements (default [])
87 * @param[in] maxiter maximum allowed number of iterations (default 100)
88 * @param[in] is_feasible true if the \p init_xs are obtained from
89 * integrating the \p init_us (rollout) (default false)
90 * @param[in] init_reg initial guess for the regularization value. Very
91 * low values are typical used with very good guess points (default 1e-9).
92 * @return A boolean that describes if convergence was reached.
93 */
94 virtual bool solve(
95 const std::vector<Eigen::VectorXd>& init_xs = DEFAULT_VECTOR,
96 const std::vector<Eigen::VectorXd>& init_us = DEFAULT_VECTOR,
97 const std::size_t maxiter = 100, const bool is_feasible = false,
98 const double reg_init = NAN) = 0;
99
100 /**
101 * @brief Compute the search direction
102 * \f$(\delta\mathbf{x}^k,\delta\mathbf{u}^k)\f$ for the current guess
103 * \f$(\mathbf{x}^k_s,\mathbf{u}^k_s)\f$.
104 *
105 * You must call `setCandidate()` first in order to define the current guess.
106 * A current guess defines a state and control trajectory
107 * \f$(\mathbf{x}^k_s,\mathbf{u}^k_s)\f$ of \f$T+1\f$ and \f$T\f$ elements,
108 * respectively.
109 *
110 * @param[in] recalc true for recalculating the derivatives at current state
111 * and control
112 * @return The search direction \f$(\delta\mathbf{x},\delta\mathbf{u})\f$ and
113 * the dual lambdas as lists of \f$T+1\f$, \f$T\f$ and \f$T+1\f$ lengths,
114 * respectively
115 */
116 virtual void computeDirection(const bool recalc) = 0;
117
118 /**
119 * @brief Try a predefined step length \f$\alpha\f$ and compute its cost
120 * improvement \f$dV\f$.
121 *
122 * It uses the search direction found by `computeDirection()` to try a
123 * determined step length \f$\alpha\f$. Therefore, it assumes that we have run
124 * `computeDirection()` first. Additionally, it returns the cost improvement
125 * \f$dV\f$ along the predefined step length \f$\alpha\f$.
126 *
127 * @param[in] steplength applied step length (\f$0\leq\alpha\leq1\f$)
128 * @return the cost improvement
129 */
130 virtual double tryStep(const double steplength = 1) = 0;
131
132 /**
133 * @brief Return a positive value that quantifies the algorithm termination
134 *
135 * These values typically represents the gradient norm which tell us that it's
136 * been reached the local minima. The stopping criteria strictly speaking
137 * depends on the search direction (calculated by `computeDirection()`) but
138 * it could also depend on the chosen step length, tested by `tryStep()`.
139 */
140 virtual double stoppingCriteria() = 0;
141
142 /**
143 * @brief Return the expected improvement \f$dV_{exp}\f$ from a given current
144 * search direction \f$(\delta\mathbf{x}^k,\delta\mathbf{u}^k)\f$
145 *
146 * For computing the expected improvement, you need to compute the search
147 * direction first via `computeDirection()`.
148 */
149 virtual const Eigen::Vector2d& expectedImprovement() = 0;
150
151 /**
152 * @brief Resizing the solver data
153 *
154 * If the shooting problem has changed after construction, then this function
155 * resizes all the data before starting resolve the problem.
156 */
157 virtual void resizeData();
158
159 /**
160 * @brief Compute the dynamic feasibility
161 * \f$\|\mathbf{f}_{\mathbf{s}}\|_{\infty,1}\f$ for the current guess
162 * \f$(\mathbf{x}^k,\mathbf{u}^k)\f$
163 *
164 * The feasibility can be computed using different norms (e.g,
165 * \f$\ell_\infty\f$ or \f$\ell_1\f$ norms). By default we use the
166 * \f$\ell_\infty\f$ norm, however, we can change the type of norm using
167 * `set_feasnorm`. Note that \f$\mathbf{f}_{\mathbf{s}}\f$ are the gaps on the
168 * dynamics, which are computed at each node as
169 * \f$\mathbf{x}^{'}-\mathbf{f}(\mathbf{x},\mathbf{u})\f$.
170 */
171 double computeDynamicFeasibility();
172
173 /**
174 * @brief Compute the feasibility of the inequality constraints for the
175 * current guess
176 *
177 * The feasibility can be computed using different norms (e.g,
178 * \f$\ell_\infty\f$ or \f$\ell_1\f$ norms). By default we use the
179 * \f$\ell_\infty\f$ norm, however, we can change the type of norm using
180 * `set_feasnorm`.
181 */
182 double computeInequalityFeasibility();
183
184 /**
185 * @brief Compute the feasibility of the equality constraints for the current
186 * guess
187 *
188 * The feasibility can be computed using different norms (e.g,
189 * \f$\ell_\infty\f$ or \f$\ell_1\f$ norms). By default we use the
190 * \f$\ell_\infty\f$ norm, however, we can change the type of norm using
191 * `set_feasnorm`.
192 */
193 double computeEqualityFeasibility();
194
195 /**
196 * @brief Set the solver candidate trajectories
197 * \f$(\mathbf{x}_s,\mathbf{u}_s)\f$
198 *
199 * The solver candidates are defined as a state and control trajectories
200 * \f$(\mathbf{x}_s,\mathbf{u}_s)\f$ of \f$T+1\f$ and \f$T\f$ elements,
201 * respectively. Additionally, we need to define the dynamic feasibility of
202 * the \f$(\mathbf{x}_s,\mathbf{u}_s)\f$ pair. Note that the trajectories are
203 * feasible if \f$\mathbf{x}_s\f$ is the resulting trajectory from the system
204 * rollout with \f$\mathbf{u}_s\f$ inputs.
205 *
206 * @param[in] xs state trajectory of \f$T+1\f$ elements (default [])
207 * @param[in] us control trajectory of \f$T\f$ elements (default [])
208 * @param[in] isFeasible true if the \p xs are obtained from integrating the
209 * \p us (rollout)
210 */
211 void setCandidate(
212 const std::vector<Eigen::VectorXd>& xs_warm = DEFAULT_VECTOR,
213 const std::vector<Eigen::VectorXd>& us_warm = DEFAULT_VECTOR,
214 const bool is_feasible = false);
215
216 /**
217 * @brief Set a list of callback functions using for the solver diagnostic
218 *
219 * Each iteration, the solver calls these set of functions in order to allowed
220 * user the diagnostic of its performance.
221 *
222 * @param callbacks set of callback functions
223 */
224 void setCallbacks(
225 const std::vector<std::shared_ptr<CallbackAbstract> >& callbacks);
226
227 /**
228 * @brief Return the list of callback functions using for diagnostic
229 */
230 const std::vector<std::shared_ptr<CallbackAbstract> >& getCallbacks() const;
231
232 /**
233 * @brief Return the shooting problem
234 */
235 const std::shared_ptr<ShootingProblem>& get_problem() const;
236
237 /**
238 * @brief Return the state trajectory \f$\mathbf{x}_s\f$
239 */
240 const std::vector<Eigen::VectorXd>& get_xs() const;
241
242 /**
243 * @brief Return the control trajectory \f$\mathbf{u}_s\f$
244 */
245 const std::vector<Eigen::VectorXd>& get_us() const;
246
247 /**
248 * @brief Return the dynamic infeasibility \f$\mathbf{f}_{s}\f$
249 */
250 const std::vector<Eigen::VectorXd>& get_fs() const;
251
252 /**
253 * @brief Return the feasibility status of the
254 * \f$(\mathbf{x}_s,\mathbf{u}_s)\f$ trajectory
255 */
256 bool get_is_feasible() const;
257
258 /**
259 * @brief Return the cost for the current guess
260 */
261 double get_cost() const;
262
263 /**
264 * @brief Return the merit for the current guess
265 */
266 double get_merit() const;
267
268 /**
269 * @brief Return the stopping-criteria value computed by `stoppingCriteria()`
270 */
271 double get_stop() const;
272
273 /**
274 * @brief Return the linear and quadratic terms of the expected improvement
275 */
276 const Eigen::Vector2d& get_d() const;
277
278 /**
279 * @brief Return the reduction in the cost function \f$\Delta V\f$
280 */
281 double get_dV() const;
282
283 /**
284 * @brief Return the reduction in the merit function \f$\Delta\Phi\f$
285 */
286 double get_dPhi() const;
287
288 /**
289 * @brief Return the expected reduction in the cost function \f$\Delta
290 * V_{exp}\f$
291 */
292 double get_dVexp() const;
293
294 /**
295 * @brief Return the expected reduction in the merit function
296 * \f$\Delta\Phi_{exp}\f$
297 */
298 double get_dPhiexp() const;
299
300 /**
301 * @brief Return the reduction in the feasibility
302 */
303 double get_dfeas() const;
304
305 /**
306 * @brief Return the total feasibility for the current guess
307 */
308 double get_feas() const;
309
310 /**
311 * @brief Return the dynamic feasibility for the current guess
312 */
313 double get_ffeas() const;
314
315 /**
316 * @brief Return the inequality feasibility for the current guess
317 */
318 double get_gfeas() const;
319
320 /**
321 * @brief Return the equality feasibility for the current guess
322 */
323 double get_hfeas() const;
324
325 /**
326 * @brief Return the dynamic feasibility for the current step length
327 */
328 double get_ffeas_try() const;
329
330 /**
331 * @brief Return the inequality feasibility for the current step length
332 */
333 double get_gfeas_try() const;
334
335 /**
336 * @brief Return the equality feasibility for the current step length
337 */
338 double get_hfeas_try() const;
339
340 /**
341 * @brief Return the primal-variable regularization
342 */
343 double get_preg() const;
344
345 /**
346 * @brief Return the dual-variable regularization
347 */
348 double get_dreg() const;
349
350 DEPRECATED("Use get_preg for primal-variable regularization",
351 double get_xreg() const;)
352 DEPRECATED("Use get_preg for primal-variable regularization",
353 double get_ureg() const;)
354
355 /**
356 * @brief Return the step length \f$\alpha\f$
357 */
358 double get_steplength() const;
359
360 /**
361 * @brief Return the threshold used for accepting a step
362 */
363 double get_th_acceptstep() const;
364
365 /**
366 * @brief Return the tolerance for stopping the algorithm
367 */
368 double get_th_stop() const;
369
370 /**
371 * @brief Return the threshold for accepting a gap as non-zero
372 */
373 double get_th_gaptol() const;
374
375 /**
376 * @brief Return the type of norm used to evaluate the dynamic and constraints
377 * feasibility
378 */
379 FeasibilityNorm get_feasnorm() const;
380
381 /**
382 * @brief Return the number of iterations performed by the solver
383 */
384 std::size_t get_iter() const;
385
386 /**
387 * @brief Modify the state trajectory \f$\mathbf{x}_s\f$
388 */
389 void set_xs(const std::vector<Eigen::VectorXd>& xs);
390
391 /**
392 * @brief Modify the control trajectory \f$\mathbf{u}_s\f$
393 */
394 void set_us(const std::vector<Eigen::VectorXd>& us);
395
396 /**
397 * @brief Modify the primal-variable regularization value
398 */
399 void set_preg(const double preg);
400
401 /**
402 * @brief Modify the dual-variable regularization value
403 */
404 void set_dreg(const double dreg);
405
406 DEPRECATED("Use set_preg for primal-variable regularization",
407 void set_xreg(const double xreg);)
408 DEPRECATED("Use set_preg for primal-variable regularization",
409 void set_ureg(const double ureg);)
410
411 /**
412 * @brief Modify the threshold used for accepting step
413 */
414 void set_th_acceptstep(const double th_acceptstep);
415
416 /**
417 * @brief Modify the tolerance for stopping the algorithm
418 */
419 void set_th_stop(const double th_stop);
420
421 /**
422 * @brief Modify the threshold for accepting a gap as non-zero
423 */
424 void set_th_gaptol(const double th_gaptol);
425
426 /**
427 * @brief Modify the current norm used for computed the dynamic and constraint
428 * feasibility
429 */
430 void set_feasnorm(const FeasibilityNorm feas_norm);
431
432 protected:
433 std::shared_ptr<ShootingProblem> problem_; //!< optimal control problem
434 std::vector<Eigen::VectorXd> xs_; //!< State trajectory
435 std::vector<Eigen::VectorXd> us_; //!< Control trajectory
436 std::vector<Eigen::VectorXd> fs_; //!< Gaps/defects between shooting nodes
437 std::vector<std::shared_ptr<CallbackAbstract> >
438 callbacks_; //!< Callback functions
439 bool is_feasible_; //!< Label that indicates is the iteration is feasible
440 bool was_feasible_; //!< Label that indicates in the previous iterate was
441 //!< feasible
442 double cost_; //!< Cost for the current guess
443 double merit_; //!< Merit for the current guess
444 double stop_; //!< Value computed by `stoppingCriteria()`
445 Eigen::Vector2d d_; //!< LQ approximation of the expected improvement
446 double dV_; //!< Reduction in the cost function computed by `tryStep()`
447 double dPhi_; //!< Reduction in the merit function computed by `tryStep()`
448 double dVexp_; //!< Expected reduction in the cost function
449 double dPhiexp_; //!< Expected reduction in the merit function
450 double dfeas_; //!< Reduction in the feasibility
451 double feas_; //!< Total feasibility for the current guess
452 double
453 ffeas_; //!< Feasibility of the dynamic constraints for the current guess
454 double gfeas_; //!< Feasibility of the inequality constraints for the current
455 //!< guess
456 double hfeas_; //!< Feasibility of the equality constraints for the current
457 //!< guess
458 double ffeas_try_; //!< Feasibility of the dynamic constraints evaluated for
459 //!< the current step length
460 double gfeas_try_; //!< Feasibility of the inequality constraints evaluated
461 //!< for the current step length
462 double hfeas_try_; //!< Feasibility of the equality constraints evaluated for
463 //!< the current step length
464 double preg_; //!< Current primal-variable regularization value
465 double dreg_; //!< Current dual-variable regularization value
466 DEPRECATED("Use preg_ for primal-variable regularization",
467 double xreg_;) //!< Current state regularization value
468 DEPRECATED("Use dreg_ for primal-variable regularization",
469 double ureg_;) //!< Current control regularization values
470 double steplength_; //!< Current applied step length
471 double th_acceptstep_; //!< Threshold used for accepting step
472 double th_stop_; //!< Tolerance for stopping the algorithm
473 double th_gaptol_; //!< Threshold limit to check non-zero gaps
474 enum FeasibilityNorm feasnorm_; //!< Type of norm used to evaluate the
475 //!< dynamics and constraints feasibility
476 std::size_t iter_; //!< Number of iteration performed by the solver
477 double tmp_feas_; //!< Temporal variables used for computed the feasibility
478 std::vector<Eigen::VectorXd> g_adj_; //!< Adjusted inequality bound
479 };
480
481 /**
482 * @brief Abstract class for solver callbacks
483 *
484 * A callback is used to diagnostic the behaviour of our solver in each
485 * iteration of it. For instance, it can be used to print values, record data or
486 * display motions.
487 */
488 class CallbackAbstract {
489 public:
490 /**
491 * @brief Initialize the callback function
492 */
493 26 CallbackAbstract() {}
494 56 virtual ~CallbackAbstract() {}
495
496 /**
497 * @brief Run the callback function given a solver
498 *
499 * @param[in] solver solver to be diagnostic
500 */
501 virtual void operator()(SolverAbstract& solver) = 0;
502 };
503
504 bool raiseIfNaN(const double value);
505
506 } // namespace crocoddyl
507
508 #endif // CROCODDYL_CORE_SOLVER_BASE_HPP_
509