Solver.hxx

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco and Enrico Bertolazzi.                                      *
 *                                                                                               *
 * The Pipal project is distributed under the MIT License.                                       *
 *                                                                                               *
 * Davide Stocco                                                               Enrico Bertolazzi *
 * University of Trento                                                     University of Trento *
 * e-mail: davide.stocco@unitn.it                             e-mail: enrico.bertolazzi@unitn.it *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#pragma once
 
#ifndef INCLUDE_PIPAL_SOLVER_HXX
#define INCLUDE_PIPAL_SOLVER_HXX
 
namespace Pipal {

  template <typename Real>
  class Solver
  {
    static_assert(std::is_floating_point_v<Real>,
      "Pipal::Solver<Real>: Real must be a floating-point type.");
 
  public:
    using ProblemPtr              = typename Problem<Real>::UniquePtr;
    using ObjectiveFunc           = typename ProblemWrapper<Real>::ObjectiveFunc;
    using ObjectiveGradientFunc   = typename ProblemWrapper<Real>::ObjectiveGradientFunc;
    using ConstraintsFunc         = typename ProblemWrapper<Real>::ConstraintsFunc;
    using ConstraintsJacobianFunc = typename ProblemWrapper<Real>::ConstraintsJacobianFunc;
    using LagrangianHessianFunc   = typename ProblemWrapper<Real>::LagrangianHessianFunc;
    using BoundsFunc              = typename ProblemWrapper<Real>::BoundsFunc;
 
  private:
    Counter          m_counter;    
    Acceptance<Real> m_acceptance; 
    Input<Real>      m_input;      
    Direction<Real>  m_direction;  
    Iterate<Real>    m_iterate;    
    Output<Real>     m_output;     
    Parameter<Real>  m_parameter;  
    ProblemPtr       m_problem;    
 
    // Some options for the solver
    bool m_verbose{false}; 
    bool m_bfgs{false};    
 
    void buildIterate();
    void evalStep();
    void updateParameters();
    void lineSearch();
    void bfgsUpdate(Vector<Real> const & s, Vector<Real> const & y);
    void updateIterate();
    void updatePoint();
    void backtracking();
    void evalFunctions();
    void evalGradients();
    void evalHessian();
    void evalNewtonMatrix();
    void evalScalings();
    void evalModels();
    void fractionToBoundary();
    void evalNewtonRhs();
    void evalSlacks();
    void evalMerit();
    void initNewtonMatrix();
    void resetDirection(Direction<Real> & d) const;
    void evalLambdaOriginal(Vector<Real> & l) const;
    Real evalViolation(Array<Real> const & cE, Array<Real> const & cI) const;
    void evalNewtonStep();
    void evalTrialStep(Direction<Real> & v) const;
    Integer checkTermination() const;
    Integer fullStepCheck();
    Integer secondOrderCorrection();
    void evalXOriginal(Vector<Real> & x);
    void setDirection(Vector<Real> const & dx, Vector<Real> const & dr1, Vector<Real> const & dr2,
      Vector<Real> const & ds1, Vector<Real> const & ds2, Vector<Real> const & dlE,
      Vector<Real> const & dlI, Real const dx_norm, Real const dl_norm);
    void evalTrialSteps(Direction<Real> & d1, Direction<Real> & d2, Direction<Real> & d3);
    void evalTrialStepCut();
    void evalLinearCombination(Direction<Real> const & d1, Direction<Real> const & d2,
      Direction<Real> const & d3, Real const a1, Real const a2, Real const a3);
    Real evalKKTError(Real const rho, Real const mu);
    void evalKKTErrors();
    void setPrimals(Vector<Real> const & x, Array<Real> const & r1, Array<Real> const & r2,
      Array<Real> const & s1, Array<Real> const & s2, Array<Real> const & lE, Array<Real> const & lI,
      Real const f, Array<Real> const & cE, Array<Real> const & cI, Real const phi);
    void buildInput(std::string const & name, Vector<Real> const & x0, Vector<Real> const & bl,
      Vector<Real> const & bu, Vector<Real> const & cl, Vector<Real> const & cu);

    void evalInfeasibility(Iterate<Real> & z) const
    {
      // Evaluate scaled and unscaled feasibility violations
      z.v  = this->evalViolation(z.cE, z.cI)/std::max(1.0, z.v0);
      z.vu = this->evalViolation(z.cEu, z.cIu);
    }

    void resetMuMaxExp() {this->m_parameter.mu_max_exp = this->m_parameter.mu_max_exp0;}

    void buildParameter(Algorithm t_algorithm) {this->m_parameter.algorithm = t_algorithm;}

    void setMuMaxExpZero() {this->m_parameter.mu_max_exp = 0.0;}

    void setRho(Real const rho) {this->m_iterate.rho = rho;}

    void setRhoLast(Real const rho) {this->m_iterate.rho_ = rho;}

    void setMu(Real const mu) {this->m_iterate.mu = mu;}

    void evalDependent()
    {
      // Evaluate quantities dependent on penalty and interior-point parameters
      this->evalSlacks();
      this->evalMerit();
      this->evalKKTErrors();
    }

    void resetCounter()
    {
      this->m_counter.f = this->m_counter.g = this->m_counter.H = this->m_counter.k = this->m_counter.M = 0;
    }

    void incrementFactorizationCount() {++this->m_counter.M;}

    void incrementFunctionCount() {++this->m_counter.f;}

    void incrementGradientCount() {++this->m_counter.g;}

    void incrementHessianCount() {++this->m_counter.H;}

    void incrementIterationCount() {++this->m_counter.k;}
 
  public:
    Solver() = default;

    Solver(
      std::string             const & name,
      ObjectiveFunc           const & objective,
      ObjectiveGradientFunc   const & objective_gradient,
      ConstraintsFunc         const & constraints,
      ConstraintsJacobianFunc const & constraints_jacobian,
      LagrangianHessianFunc   const & lagrangian_hessian,
      BoundsFunc              const & primal_lower_bounds,
      BoundsFunc              const & primal_upper_bounds,
      BoundsFunc              const & constraints_lower_bounds,
      BoundsFunc              const & constraints_upper_bounds
    ) : m_problem(std::make_unique<ProblemWrapper<Real>>(
      name,
      objective,
      objective_gradient,
      constraints,
      constraints_jacobian,
      lagrangian_hessian,
      primal_lower_bounds,
      primal_upper_bounds,
      constraints_lower_bounds,
      constraints_upper_bounds
    )) {}

    Solver(ProblemPtr && problem) : m_problem(std::move(problem)) {}

    Solver(Solver const &) = delete;

    Solver & operator=(Solver const &) = delete;

    Solver(Solver &&) = delete;

    Solver & operator=(Solver &&) = delete;

    ~Solver() = default;

    void problem(ProblemPtr && problem) {this->m_problem = std::move(problem);}

    Problem<Real> const & problem() const {return this->m_problem.get();}

    bool verbose_mode() const {return this->m_verbose;}

    void verbose_mode(bool const t_verbose) {this->m_verbose = t_verbose;}

    bool bfgs() const {return this->m_bfgs;}

    void bfgs(bool const t_bfgs) {this->m_bfgs = t_bfgs;}

    Algorithm algorithm() const {return this->m_parameter.algorithm;}

    void algorithm(Algorithm const t_algorithm) {this->m_parameter.algorithm = t_algorithm;}

    void tolerance(Real const t_tolerance)
    {
      PIPAL_ASSERT(t_tolerance > 0.0,
        "Pipal::Solver::tolerance(...): input value must be positive");
      this->m_parameter.opt_err_tol = t_tolerance;
    }

    Real tolerance() const {return this->m_parameter.opt_err_tol;}

    void max_iterations(Integer const t_max_iterations) {
      PIPAL_ASSERT(t_max_iterations > 0,
        "Pipal::Solver::max_iterations(...): input value must be positive");
      this->m_parameter.iter_max = t_max_iterations;
    }

    Integer max_iterations() const {return this->m_parameter.iter_max;}

    bool optimize(Vector<Real> const & x_guess, Vector<Real> & x_sol)
    {
 
      #define CMD "Pipal::Solver::optimize(...): "
 
      // Create alias for easier access
      Counter          & c{this->m_counter};
      Input<Real>      & i{this->m_input};
      Iterate<Real>    & z{this->m_iterate};
      Direction<Real>  & d{this->m_direction};
      Acceptance<Real> & a{this->m_acceptance};
 
      // Check that the problem is set
      PIPAL_ASSERT(this->m_problem.get() != nullptr,
        CMD "problem not set, use 'problem(...)' method to set it");
 
      // Get variable bounds
      Vector<Real> bl, bu;
      PIPAL_ASSERT(this->m_problem->primal_lower_bounds(bl),
        CMD "error in evaluating lower bounds on primal variables");
      PIPAL_ASSERT(this->m_problem->primal_upper_bounds(bu),
        CMD "error in evaluating upper bounds on primal variables");
 
      // Get constraint bounds
      Vector<Real> cl, cu;
      PIPAL_ASSERT(this->m_problem->constraints_lower_bounds(cl),
        CMD "error in evaluating lower bounds on constraints");
      PIPAL_ASSERT(this->m_problem->constraints_upper_bounds(cu),
        CMD "error in evaluating upper bounds on constraints");
 
      // Reset counters
      resetCounter();
 
      // Fill input structure
      buildInput(this->m_problem->name(), x_guess, bl, bu, cl, cu);
      buildIterate();
      resetDirection(d);
 
      // Print header and break line
      if (this->m_verbose) {this->m_output.printHeader(i, z); this->m_output.printBreak(c);}
 
      // Iterations loop
      while (!this->checkTermination()) {
 
        // Print iterate
        if (this->m_verbose) {this->m_output.printIterate(c, z);}
 
        // Evaluate the step
        this->evalStep();
 
        // Print direction
        if (this->m_verbose) {this->m_output.printDirection(z, d);}
 
        this->lineSearch();
 
        // Print accepted
        if (this->m_verbose) {this->m_output.printAcceptance(a);}
 
        this->updateIterate();
 
        // Increment iteration counter
        this->incrementIterationCount();
 
        // Print break
        if (this->m_verbose) {this->m_output.printBreak(c);}
      }
      // Print footer and terminate
      if (this->m_verbose) {this->m_output.printFooter(c, z, this->checkTermination());}
 
      // Get solution in original variables
      this->evalXOriginal(x_sol);
 
      // Return success if is finite
      return x_sol.allFinite();
 
      #undef CMD
    }

    void getSolution(Vector<Real> & x, Vector<Real> & l)
    {
      this->evalXOriginal(x);
      this->evalLambdaOriginal(l);
    }
 
  }; // class Solver
 
} // namespace Solver
 
#endif // INCLUDE_PIPAL_SOLVER_HXX