qpmad
Eigen-based C++ QP solver.
Loading...
Searching...
No Matches
solver.h
Go to the documentation of this file.
1/**
2 @file
3 @author Alexander Sherikov
4
5 @copyright 2017 Alexander Sherikov. Licensed under the Apache License,
6 Version 2.0. (see LICENSE or http://www.apache.org/licenses/LICENSE-2.0)
7
8 @brief
9*/
10
11
12#pragma once
13
14#include "common.h"
15#include "givens.h"
16#include "input_parser.h"
17#include "inverse.h"
18#include "solver_parameters.h"
19#include "constraint_status.h"
20#include "chosen_constraint.h"
21#include "active_set.h"
22#include "factorization_data.h"
23
24#ifdef QPMAD_ENABLE_TRACING
25# include "testing.h"
26#endif
27
28// is_eigen_type requires Boost
29#define QPMAD_UTILS_EIGEN_VECTOR_ENABLER(Type) const typename std::enable_if<Type::IsVectorAtCompileTime>::type * = NULL
30#define QPMAD_UTILS_EIGEN_MATRIX_ENABLER(Type) \
31 const typename std::enable_if<(Type::ColsAtCompileTime == Eigen::Dynamic or Type::ColsAtCompileTime > 1)>::type \
32 * = NULL
33
34
35namespace qpmad
36{
37 template <typename t_Scalar, int t_primal_size, int t_has_bounds, int t_general_ctr_number>
39 {
40 public:
42 {
43 OK = 0,
45 };
46
47 template <int t_rows>
48 using Vector = Eigen::Matrix<t_Scalar, t_rows, 1>;
49 template <int t_rows, int t_cols>
50 using Matrix = Eigen::Matrix<t_Scalar, t_rows, t_cols>; /// @deprecated
51 using Scalar = t_Scalar;
52
53
54 protected:
55 const struct InputPlaceholders
56 {
57 Eigen::Matrix<t_Scalar, Eigen::Dynamic, 1> empty_vector_;
58 Eigen::Matrix<t_Scalar, Eigen::Dynamic, Eigen::Dynamic> empty_matrix_;
61
62
63 protected:
65 Eigen::Dynamic == t_general_ctr_number ?
66 Eigen::Dynamic :
67 (0 == t_has_bounds ? t_general_ctr_number :
68 (Eigen::Dynamic == t_primal_size ? Eigen::Dynamic :
69 t_general_ctr_number + t_primal_size));
72
75
77
80
82
83 Eigen::Array<uint8_t, num_constraints_compile_time_, 1> constraints_status_;
84
86
87 std::ptrdiff_t iter_counter_;
88
90
91
92 public:
94 {
95 iter_counter_ = 0;
98 }
99
100
101 /**
102 * @brief Returns type of the Hessian produced by the latest execution
103 * of `solve()`.
104 */
106 {
107 return (hessian_type_);
108 }
109
110
111 /**
112 * @brief Returns number of inequality iterations during the latest
113 * execution of `solve()`.
114 */
115 std::ptrdiff_t getNumberOfInequalityIterations() const
116 {
117 return (iter_counter_);
118 }
119
120
121 /**
122 * @brief Returns dual variables (Lagrange multipliers) corresponding
123 * to inequality constraints. Must be called after successful
124 * `solve()`, the result is undefined if previous call to `solve()`
125 * failed.
126 *
127 * @tparam t_status_size
128 * @tparam t_dual_size
129 * @tparam t_index_size
130 *
131 * @param[out] dual dual variables
132 * @param[out] indices constraint indices corresponding to the dual
133 * variables, index 0 corresponds to the first simple bound if present
134 * or to the first general constraint otherwise
135 * @param[out] is_lower flags indicating if lower or upper bound is
136 * active
137 */
138 template <int t_status_size, int t_dual_size, int t_index_size>
141 Eigen::Matrix<MatrixIndex, t_index_size, 1> &indices,
142 Eigen::Matrix<bool, t_status_size, 1> &is_lower) const
143 {
145 {
146 const MatrixIndex size = active_set_.size_ - active_set_.num_equalities_;
147
148 dual.resize(size);
149 indices.resize(size);
150 is_lower.resize(size);
151
152 for (MatrixIndex i = active_set_.num_equalities_; i < active_set_.size_; ++i)
153 {
154 const std::size_t output_index = i - active_set_.num_equalities_;
155
156 dual(output_index) = dual_(i);
157 indices(output_index) = active_set_.getIndex(i);
158 is_lower(output_index) =
160 }
161 }
162 else
163 {
164 dual.resize(0);
165 indices.resize(0);
166 is_lower.resize(0);
167 }
168 }
169
170
171 template <class t_primal, class t_H, class t_h, class t_A, class t_Alb, class t_Aub>
173 t_primal &primal,
174 t_H &H,
175 const t_h &h,
176 const t_A &A,
177 const t_Alb &Alb,
178 const t_Aub &Aub,
179 // accept only Eigen types as inputs
186 {
187 return (
188 solve(primal,
189 H,
190 h,
193 A,
194 Alb,
195 Aub,
197 }
198
199
200 template <class t_primal, class t_H, class t_h, class t_A, class t_Alb>
202 t_primal &primal,
203 t_H &H,
204 const t_h &h,
205 const t_A &A,
206 const t_Alb &Alb,
207 // accept only Eigen types as inputs
213 {
214 return (
215 solve(primal,
216 H,
217 h,
220 A,
221 Alb,
222 Alb,
224 }
225
226
227 template <class t_primal, class t_H, class t_h, class t_lb, class t_ub, class t_A, class t_Alb, class t_Aub>
229 t_primal &primal,
230 t_H &H,
231 const t_h &h,
232 const t_lb &lb,
233 const t_ub &ub,
234 const t_A &A,
235 const t_Alb &Alb,
236 const t_Aub &Aub,
237 // accept only Eigen types as inputs
246 {
247 return (solve(primal, H, h, lb, ub, A, Alb, Aub, input_placeholders_.solver_parameters_));
248 }
249
250
251 template <class t_primal, class t_H, class t_h, class t_lb, class t_ub, class t_A, class t_Alb>
253 t_primal &primal,
254 t_H &H,
255 const t_h &h,
256 const t_lb &lb,
257 const t_ub &ub,
258 const t_A &A,
259 const t_Alb &Alb,
260 // accept only Eigen types as inputs
268 {
269 return (solve(primal, H, h, lb, ub, A, Alb, Alb, input_placeholders_.solver_parameters_));
270 }
271
272
273 template <class t_primal, class t_H, class t_h, class t_lb, class t_ub>
275 t_primal &primal,
276 t_H &H,
277 const t_h &h,
278 const t_lb &lb,
279 const t_ub &ub,
280 const SolverParameters &param,
281 // accept only Eigen types as inputs
287 {
288 return (
289 solve(primal,
290 H,
291 h,
292 lb,
293 ub,
297 param));
298 }
299
300
301 template <class t_primal, class t_H, class t_h, class t_lb, class t_ub>
303 t_primal &primal,
304 t_H &H,
305 const t_h &h,
306 const t_lb &lb,
307 const t_ub &ub,
308 // accept only Eigen types as inputs
314 {
315 return (
316 solve(primal,
317 H,
318 h,
319 lb,
320 ub,
325 }
326
327
328
329 template <class t_primal, class t_H, class t_h, class t_lb, class t_ub, class t_A, class t_Alb, class t_Aub>
331 t_primal &primal,
332 t_H &H,
333 const t_h &h,
334 const t_lb &lb,
335 const t_ub &ub,
336 const t_A &A,
337 const t_Alb &Alb,
338 const t_Aub &Aub,
339 const SolverParameters &param,
340 // accept only Eigen types as inputs
349 {
350 QPMAD_TRACE(std::setprecision(std::numeric_limits<double>::digits10));
351
353 iter_counter_ = 0;
354
355 parseObjective(H, h);
356 parseSimpleBounds(lb, ub);
357 parseGeneralConstraints(A, Alb, Aub);
358
359
361
362 if (0 == h_size_)
363 {
364 // trivial unconstrained optimum
365 primal.setZero(primal_size_);
366
367 if (0 == num_constraints_)
368 {
369 // trivial solution
370 return (OK);
371 }
372 }
373
374
375 switch (param.hessian_type_)
376 {
379 /* Falls through. */
382 // unconstrained optimum
383 if (h_size_ > 0)
384 {
385 primal = H.template triangularView<Eigen::Lower>().solve(-h);
386 H.transpose().template triangularView<Eigen::Upper>().solveInPlace(primal);
387 }
388 break;
389
392 // unconstrained optimum
393 if (h_size_ > 0)
394 {
395 primal_step_direction_.noalias() = H.template triangularView<Eigen::Upper>().transpose() * -h;
396 primal.noalias() = H.template triangularView<Eigen::Upper>() * primal_step_direction_;
397 }
398 break;
399
400 default:
401 QPMAD_UTILS_THROW("Malformed solver parameters!");
402 break;
403 }
404
405
406 if (0 == num_constraints_)
407 {
408 // return unconstrained optimum
409 return (OK);
410 }
411
412
413 // check consistency of general constraints and activate
414 // equality constraints
416 MatrixIndex num_equalities = 0;
417 for (MatrixIndex i = 0; i < num_constraints_; ++i)
418 {
420
421 double lb_i;
422 double ub_i;
423
424 if (true == chosen_ctr_.is_simple_)
425 {
426 lb_i = lb(i);
427 ub_i = ub(i);
428 }
429 else
430 {
434 }
435
436
437 if (lb_i - param.tolerance_ > ub_i)
438 {
440 QPMAD_UTILS_THROW("Inconsistent constraints (lb > ub).");
441 }
442
443 if (std::abs(lb_i - ub_i) > param.tolerance_)
444 {
446 }
447 else
448 {
450 ++num_equalities;
451
452
453 if (true == chosen_ctr_.is_simple_)
454 {
455 chosen_ctr_.violation_ = lb_i - primal(i);
456 }
457 else
458 {
460 }
461
463
464
465 // if 'primal_size_' constraints are already activated
466 // all other constraints are linearly dependent
467 if (active_set_.hasEmptySpace())
468 {
469 double ctr_i_dot_primal_step_direction;
470
471 if (true == chosen_ctr_.is_simple_)
472 {
473 factorization_data_.computeEqualityPrimalStep(primal_step_direction_, i, active_set_.size_);
474
475 ctr_i_dot_primal_step_direction = primal_step_direction_(i);
476 }
477 else
478 {
479 factorization_data_.computeEqualityPrimalStep(
482 active_set_.size_);
483
484 ctr_i_dot_primal_step_direction =
486 }
487
488 // if step direction is a zero vector, constraint is
489 // linearly dependent with previously added constraints
490 if (ctr_i_dot_primal_step_direction < -param.tolerance_)
491 {
492 double primal_step_length_ = chosen_ctr_.violation_ / ctr_i_dot_primal_step_direction;
493
494 primal.noalias() += primal_step_length_ * primal_step_direction_;
495
496 if (false
497 == factorization_data_.update(
499 {
500 QPMAD_UTILS_THROW("Failed to add an equality constraint -- is this possible?");
501 }
502 active_set_.addEquality(i);
503
504 continue;
505 }
506 } // otherwise -- linear dependence
507
508 // this point is reached if constraint is linearly dependent
509
510 // check if this constraint is actually satisfied
511 if (std::abs(chosen_ctr_.violation_) > param.tolerance_)
512 {
513 // nope it is not
515 QPMAD_UTILS_THROW("Infeasible equality constraints");
516 }
517 // otherwise keep going
518 }
519 }
520
521
522 if (num_equalities == num_constraints_)
523 {
524 // exit early -- avoid unnecessary memory allocations
525 return (OK);
526 }
527
528
529 ReturnStatus return_status;
530 chooseConstraint(primal, lb, ub, A, Alb, Aub, param.tolerance_);
531
532
533 if (std::abs(chosen_ctr_.violation_) < param.tolerance_)
534 {
535 // all constraints are satisfied
536 return_status = OK;
537 }
538 else
539 {
540 return_status = MAXIMAL_NUMBER_OF_ITERATIONS;
541
542
545
546
547 double chosen_ctr_dot_primal_step_direction = 0.0;
548
549 //
550 factorization_data_.computeInequalityDualStep(dual_step_direction_, chosen_ctr_, A, active_set_);
551 if (active_set_.hasEmptySpace())
552 {
553 // compute step direction in primal space
554 factorization_data_.computeInequalityPrimalStep(primal_step_direction_, active_set_);
555 chosen_ctr_dot_primal_step_direction =
557 }
558
559
560 // last iteration is not counted, so iter_counter_ starts with 1.
561 for (iter_counter_ = 1; (param.max_iter_ < 0) or (iter_counter_ <= param.max_iter_); ++iter_counter_)
562 {
563 QPMAD_TRACE(">>>>>>>>>" << iter_counter_ << "<<<<<<<<<");
564#ifdef QPMAD_ENABLE_TRACING
565 testing::computeObjective(H, h, primal);
566#endif
567 QPMAD_TRACE("||| Chosen ctr index = " << chosen_ctr_.index_);
568 QPMAD_TRACE("||| Chosen ctr dual = " << chosen_ctr_.dual_);
569 QPMAD_TRACE("||| Chosen ctr violation = " << chosen_ctr_.violation_);
570
571
572 // check dual feasibility
573 MatrixIndex dual_blocking_index = primal_size_;
574 double dual_step_length = std::numeric_limits<double>::infinity();
575 for (MatrixIndex i = active_set_.num_equalities_; i < active_set_.size_; ++i)
576 {
577 if (dual_step_direction_(i) < -param.tolerance_)
578 {
579 const double dual_step_length_i = -dual_(i) / dual_step_direction_(i);
580 if (dual_step_length_i < dual_step_length)
581 {
582 dual_step_length = dual_step_length_i;
583 dual_blocking_index = i;
584 }
585 }
586 }
587
588
589#ifdef QPMAD_ENABLE_TRACING
591 H,
592 h,
593 primal,
594 A,
598 dual_,
600#endif
601
602
603 if (active_set_.hasEmptySpace()
604 // if step direction is a zero vector, constraint is
605 // linearly dependent with previously added constraints
606 && (std::abs(chosen_ctr_dot_primal_step_direction) > param.tolerance_))
607 {
608 double step_length = -chosen_ctr_.violation_ / chosen_ctr_dot_primal_step_direction;
609
610 QPMAD_TRACE("======================");
611 QPMAD_TRACE("||| Primal step length = " << step_length);
612 QPMAD_TRACE("||| Dual step length = " << dual_step_length);
613 QPMAD_TRACE("======================");
614
615
616 bool partial_step = false;
618 (step_length >= 0.0) && (dual_step_length >= 0.0),
619 "Non-negative step lengths expected.");
620 if (dual_step_length <= step_length)
621 {
622 step_length = dual_step_length;
623 partial_step = true;
624 }
625
626
627 primal.noalias() += step_length * primal_step_direction_;
628
629 dual_.segment(active_set_.num_equalities_, active_set_.num_inequalities_).noalias() +=
630 step_length
631 * dual_step_direction_.segment(
632 active_set_.num_equalities_, active_set_.num_inequalities_);
633 chosen_ctr_.dual_ += step_length;
634 chosen_ctr_.violation_ += step_length * chosen_ctr_dot_primal_step_direction;
635
636 QPMAD_TRACE("||| Chosen ctr dual = " << chosen_ctr_.dual_);
637 QPMAD_TRACE("||| Chosen ctr violation = " << chosen_ctr_.violation_);
638
639
640 if ((partial_step)
641 // if violation is almost zero -- assume that a full step is made
642 && (std::abs(chosen_ctr_.violation_) > param.tolerance_))
643 {
644 QPMAD_TRACE("||| PARTIAL STEP");
645 // deactivate blocking constraint
646 constraints_status_[active_set_.getIndex(dual_blocking_index)] = ConstraintStatus::INACTIVE;
647
648 dropElementWithoutResize(dual_, dual_blocking_index, active_set_.size_);
649
650 factorization_data_.downdate(dual_blocking_index, active_set_.size_);
651
652 active_set_.removeInequality(dual_blocking_index);
653
654 // compute step direction in primal & dual space
655 factorization_data_.updateStepsAfterPartialStep(
657 chosen_ctr_dot_primal_step_direction =
659 }
660 else
661 {
662 QPMAD_TRACE("||| FULL STEP");
663 // activate constraint
664 if (false
665 == factorization_data_.update(
667 {
668 QPMAD_UTILS_THROW("Failed to add an inequality constraint -- is this possible?");
669 }
670
672 {
674 }
675 else
676 {
678 }
680 active_set_.addInequality(chosen_ctr_.index_);
681
682 chooseConstraint(primal, lb, ub, A, Alb, Aub, param.tolerance_);
683
684 if (std::abs(chosen_ctr_.violation_) < param.tolerance_)
685 {
686 // all constraints are satisfied
687 return_status = OK;
688 break;
689 }
690
691 chosen_ctr_dot_primal_step_direction = 0.0;
692 factorization_data_.computeInequalityDualStep(
694 if (active_set_.hasEmptySpace())
695 {
696 // compute step direction in primal & dual space
697 factorization_data_.computeInequalityPrimalStep(primal_step_direction_, active_set_);
698 chosen_ctr_dot_primal_step_direction =
700 }
701 }
702 }
703 else
704 {
705 if (dual_blocking_index == primal_size_)
706 {
707 QPMAD_UTILS_THROW("Infeasible inequality constraints.");
708 }
709 else
710 {
711 QPMAD_TRACE("======================");
712 QPMAD_TRACE("||| Dual step length = " << dual_step_length);
713 QPMAD_TRACE("======================");
714
715 // otherwise -- deactivate
716 dual_.segment(active_set_.num_equalities_, active_set_.num_inequalities_).noalias() +=
717 dual_step_length
718 * dual_step_direction_.segment(
719 active_set_.num_equalities_, active_set_.num_inequalities_);
720 chosen_ctr_.dual_ += dual_step_length;
721
722 constraints_status_[active_set_.getIndex(dual_blocking_index)] = ConstraintStatus::INACTIVE;
723
724 dropElementWithoutResize(dual_, dual_blocking_index, active_set_.size_);
725
726 factorization_data_.downdate(dual_blocking_index, active_set_.size_);
727
728 active_set_.removeInequality(dual_blocking_index);
729
730 // compute step direction in primal & dual space
731 factorization_data_.updateStepsAfterPureDualStep(
733 chosen_ctr_dot_primal_step_direction =
735 }
736 }
737 }
738 }
739
740#ifdef QPMAD_ENABLE_TRACING
742 {
744
747 }
748 else
749 {
750 QPMAD_TRACE("||| NO ACTIVE CONSTRAINTS");
751 }
752#endif
753 return (return_status);
754 }
755
756
758 const MatrixIndex primal_size,
759 const MatrixIndex num_simple_bounds,
760 const MatrixIndex num_general_constraints)
761 {
762 reserveMachinery(primal_size, num_general_constraints);
763 constraints_status_.resize(num_simple_bounds + num_general_constraints);
764 reserveDual(primal_size);
765 }
766
767
768 private:
769 void reserveMachinery(const MatrixIndex primal_size, const MatrixIndex num_general_constraints)
770 {
771 active_set_.initialize(primal_size);
772 primal_step_direction_.resize(primal_size);
773 general_ctr_dot_primal_.resize(num_general_constraints);
774
775 factorization_data_.reserve(primal_size);
776 }
777
778
779 void reserveDual(const MatrixIndex primal_size)
780 {
781 dual_.resize(primal_size);
782 dual_step_direction_.resize(primal_size);
783 }
784
785
786 template <class t_MatrixType>
787 void initializeMachineryLazy(t_MatrixType &H, const bool return_inverted_cholesky_factor)
788 {
790 {
792
793 factorization_data_.initialize(H, hessian_type_, primal_size_, return_inverted_cholesky_factor);
794 if (return_inverted_cholesky_factor)
795 {
797 }
798
800 }
801 }
802
803
804 template <
805 class t_Primal,
806 class t_LowerBounds,
807 class t_UpperBounds,
808 class t_Constraints,
809 class t_ConstraintsLowerBounds,
810 class t_ConstraintsUpperBounds>
812 const t_Primal &primal,
813 const t_LowerBounds &lb,
814 const t_UpperBounds &ub,
815 const t_Constraints &A,
816 const t_ConstraintsLowerBounds &Alb,
817 const t_ConstraintsUpperBounds &Aub,
818 const double tolerance)
819 {
821
822
823 for (MatrixIndex i = 0; i < num_simple_bounds_; ++i)
824 {
826 {
827 checkConstraintViolation(i, lb(i), ub(i), primal(i));
828 }
829 }
830
831 if ((std::abs(chosen_ctr_.violation_) < tolerance) && (num_general_constraints_ > 0))
832 {
833 general_ctr_dot_primal_.noalias() = A * primal;
835 {
837 {
839 i,
840 Alb(i - num_simple_bounds_),
841 Aub(i - num_simple_bounds_),
843 }
844 }
846 {
848 }
849 }
850
853 }
854
855
857 const MatrixIndex i,
858 const double lb_i,
859 const double ub_i,
860 const double ctr_i_dot_primal)
861 {
862 double ctr_violation_i = ctr_i_dot_primal - lb_i;
863 if (ctr_violation_i < -std::abs(chosen_ctr_.violation_))
864 {
865 chosen_ctr_.violation_ = ctr_violation_i;
867 }
868 else
869 {
870 ctr_violation_i = ctr_i_dot_primal - ub_i;
871 if (ctr_violation_i > std::abs(chosen_ctr_.violation_))
872 {
873 chosen_ctr_.violation_ = ctr_violation_i;
875 }
876 }
877 }
878
879
880 template <class t_VectorType, class t_MatrixType>
881 double getConstraintDotPrimalStepDirection(const t_VectorType &primal_step_direction, const t_MatrixType &A)
882 const
883 {
885 {
886 return (primal_step_direction(chosen_ctr_.index_));
887 }
888 else
889 {
890 return (A.row(chosen_ctr_.general_constraint_index_) * primal_step_direction);
891 }
892 }
893
894
895 template <int t_rows, int t_cols>
896 static void factorizeCholeskyInPlace(Eigen::Matrix<t_Scalar, t_rows, t_cols> &H)
897 {
898 const Eigen::LLT<Eigen::Ref<typename std::decay<decltype(H)>::type>, Eigen::Lower> llt(H);
900 Eigen::Success == llt.info(), "Could not perform Cholesky decomposition of the Hessian (dense).");
901 }
902
903
904 template <int t_Options, typename t_StorageIndex>
905 static void factorizeCholeskyInPlace(Eigen::SparseMatrix<t_Scalar, t_Options, t_StorageIndex> &H)
906 {
907 const Eigen::SimplicialLLT<typename std::decay<decltype(H)>::type, Eigen::Lower> llt(H);
909 Eigen::Success == llt.info(), "Could not perform Cholesky decomposition of the Hessian (sparse).");
910 H = llt.matrixL();
911 }
912 };
913
914
916} // namespace qpmad
MatrixIndex general_constraint_index_
MatrixIndex num_simple_bounds_
Definition: input_parser.h:20
void parseObjective(const t_MatrixType &H, const t_VectorType &h)
Definition: input_parser.h:35
MatrixIndex h_size_
Definition: input_parser.h:19
MatrixIndex num_general_constraints_
Definition: input_parser.h:21
void parseSimpleBounds(const t_VectorTypelb &lb, const t_VectorTypeub &ub)
Definition: input_parser.h:48
void parseGeneralConstraints(const t_MatrixTypeA &A, const t_VectorTypelb &lb, const t_VectorTypeub &ub)
Definition: input_parser.h:68
MatrixIndex primal_size_
Definition: input_parser.h:18
static void factorizeCholeskyInPlace(Eigen::Matrix< t_Scalar, t_rows, t_cols > &H)
Definition: solver.h:896
static const MatrixIndex num_constraints_compile_time_
Definition: solver.h:64
void initializeMachineryLazy(t_MatrixType &H, const bool return_inverted_cholesky_factor)
Definition: solver.h:787
void reserveDual(const MatrixIndex primal_size)
Definition: solver.h:779
ReturnStatus solve(t_primal &primal, t_H &H, const t_h &h, const t_lb &lb, const t_ub &ub, const t_A &A, const t_Alb &Alb, QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_primal), QPMAD_UTILS_EIGEN_MATRIX_ENABLER(t_H), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_h), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_lb), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_ub), QPMAD_UTILS_EIGEN_MATRIX_ENABLER(t_A), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_Alb))
Definition: solver.h:252
const struct qpmad::SolverTemplate::InputPlaceholders input_placeholders_
FactorizationData< t_Scalar, t_primal_size > factorization_data_
Definition: solver.h:74
Vector< t_primal_size > dual_
Definition: solver.h:76
void checkConstraintViolation(const MatrixIndex i, const double lb_i, const double ub_i, const double ctr_i_dot_primal)
Definition: solver.h:856
SolverParameters::HessianType getHessianType() const
Returns type of the Hessian produced by the latest execution of solve().
Definition: solver.h:105
ReturnStatus solve(t_primal &primal, t_H &H, const t_h &h, const t_lb &lb, const t_ub &ub, QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_primal), QPMAD_UTILS_EIGEN_MATRIX_ENABLER(t_H), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_h), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_lb), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_ub))
Definition: solver.h:302
SolverParameters::HessianType hessian_type_
Definition: solver.h:89
ReturnStatus solve(t_primal &primal, t_H &H, const t_h &h, const t_lb &lb, const t_ub &ub, const t_A &A, const t_Alb &Alb, const t_Aub &Aub, QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_primal), QPMAD_UTILS_EIGEN_MATRIX_ENABLER(t_H), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_h), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_lb), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_ub), QPMAD_UTILS_EIGEN_MATRIX_ENABLER(t_A), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_Alb), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_Aub))
Definition: solver.h:228
Vector< t_primal_size > dual_step_direction_
Definition: solver.h:79
static void factorizeCholeskyInPlace(Eigen::SparseMatrix< t_Scalar, t_Options, t_StorageIndex > &H)
Definition: solver.h:905
ReturnStatus solve(t_primal &primal, t_H &H, const t_h &h, const t_A &A, const t_Alb &Alb, const t_Aub &Aub, QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_primal), QPMAD_UTILS_EIGEN_MATRIX_ENABLER(t_H), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_h), QPMAD_UTILS_EIGEN_MATRIX_ENABLER(t_A), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_Alb), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_Aub))
Definition: solver.h:172
@ MAXIMAL_NUMBER_OF_ITERATIONS
Definition: solver.h:44
ReturnStatus solve(t_primal &primal, t_H &H, const t_h &h, const t_A &A, const t_Alb &Alb, QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_primal), QPMAD_UTILS_EIGEN_MATRIX_ENABLER(t_H), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_h), QPMAD_UTILS_EIGEN_MATRIX_ENABLER(t_A), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_Alb))
Definition: solver.h:201
void chooseConstraint(const t_Primal &primal, const t_LowerBounds &lb, const t_UpperBounds &ub, const t_Constraints &A, const t_ConstraintsLowerBounds &Alb, const t_ConstraintsUpperBounds &Aub, const double tolerance)
Definition: solver.h:811
t_Scalar Scalar
Definition: solver.h:51
Eigen::Matrix< t_Scalar, t_rows, 1 > Vector
Definition: solver.h:48
std::ptrdiff_t iter_counter_
Definition: solver.h:87
void reserveMachinery(const MatrixIndex primal_size, const MatrixIndex num_general_constraints)
Definition: solver.h:769
MatrixIndex num_constraints_
Definition: solver.h:70
ChosenConstraint chosen_ctr_
Definition: solver.h:85
std::ptrdiff_t getNumberOfInequalityIterations() const
Returns number of inequality iterations during the latest execution of solve().
Definition: solver.h:115
Vector< t_general_ctr_number > general_ctr_dot_primal_
Definition: solver.h:81
ActiveSet< t_primal_size > active_set_
Definition: solver.h:73
void reserve(const MatrixIndex primal_size, const MatrixIndex num_simple_bounds, const MatrixIndex num_general_constraints)
Definition: solver.h:757
Eigen::Array< uint8_t, num_constraints_compile_time_, 1 > constraints_status_
Definition: solver.h:83
ReturnStatus solve(t_primal &primal, t_H &H, const t_h &h, const t_lb &lb, const t_ub &ub, const t_A &A, const t_Alb &Alb, const t_Aub &Aub, const SolverParameters &param, QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_primal), QPMAD_UTILS_EIGEN_MATRIX_ENABLER(t_H), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_h), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_lb), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_ub), QPMAD_UTILS_EIGEN_MATRIX_ENABLER(t_A), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_Alb), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_Aub))
Definition: solver.h:330
Vector< t_primal_size > primal_step_direction_
Definition: solver.h:78
bool machinery_initialized_
Definition: solver.h:71
ReturnStatus solve(t_primal &primal, t_H &H, const t_h &h, const t_lb &lb, const t_ub &ub, const SolverParameters &param, QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_primal), QPMAD_UTILS_EIGEN_MATRIX_ENABLER(t_H), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_h), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_lb), QPMAD_UTILS_EIGEN_VECTOR_ENABLER(t_ub))
Definition: solver.h:274
void getInequalityDual(Vector< t_dual_size > &dual, Eigen::Matrix< MatrixIndex, t_index_size, 1 > &indices, Eigen::Matrix< bool, t_status_size, 1 > &is_lower) const
Returns dual variables (Lagrange multipliers) corresponding to inequality constraints....
Definition: solver.h:139
double getConstraintDotPrimalStepDirection(const t_VectorType &primal_step_direction, const t_MatrixType &A) const
Definition: solver.h:881
Eigen::Matrix< t_Scalar, t_rows, t_cols > Matrix
Definition: solver.h:50
#define QPMAD_TRACE(info)
Definition: common.h:27
#define QPMAD_UTILS_PERSISTENT_ASSERT(condition, message)
#define QPMAD_UTILS_THROW(s)
#define QPMAD_UTILS_ASSERT(condition, message)
double computeObjective(const Eigen::MatrixXd &H, const Eigen::VectorXd &h, const Eigen::VectorXd &primal)
Definition: testing.h:21
void printActiveSet(const t_ActiveSet &active_set, const t_ConstraintStatuses &constraints_status, const Eigen::VectorXd &dual)
Definition: testing.h:165
void checkLagrangeMultipliers(const Eigen::MatrixXd &H, const Eigen::VectorXd &h, const Eigen::VectorXd &primal, const Eigen::MatrixXd &A, const t_ActiveSet &active_set, const MatrixIndex &num_simple_bounds, const t_ConstraintStatuses &constraints_status, const Eigen::VectorXd &dual, const Eigen::VectorXd &dual_direction=Eigen::VectorXd())
Definition: testing.h:39
void dropElementWithoutResize(t_VectorType &vector, const MatrixIndex index, const MatrixIndex size)
Definition: common.h:40
EIGEN_DEFAULT_DENSE_INDEX_TYPE MatrixIndex
Definition: common.h:33
#define QPMAD_UTILS_EIGEN_MATRIX_ENABLER(Type)
Definition: solver.h:30
#define QPMAD_UTILS_EIGEN_VECTOR_ENABLER(Type)
Definition: solver.h:29
Eigen::Matrix< t_Scalar, Eigen::Dynamic, 1 > empty_vector_
Definition: solver.h:57
Eigen::Matrix< t_Scalar, Eigen::Dynamic, Eigen::Dynamic > empty_matrix_
Definition: solver.h:58