GCC Code Coverage Report


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