ipopt-src 0.2.3+3.14.16

// Copyright (C) 2004, 2006 International Business Machines and others.
// All Rights Reserved.
// This code is published under the Eclipse Public License.
//
// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13

#ifndef __IPRESTOIPOPTNLP_HPP__
#define __IPRESTOIPOPTNLP_HPP__

#include "IpIpoptNLP.hpp"
#include "IpIpoptData.hpp"
#include "IpIpoptCalculatedQuantities.hpp"
#include "IpCompoundMatrix.hpp"
#include "IpCompoundSymMatrix.hpp"
#include "IpCompoundVector.hpp"
#include "IpIdentityMatrix.hpp"
#include "IpDiagMatrix.hpp"
#include "IpZeroMatrix.hpp"
#include "IpOrigIpoptNLP.hpp"

namespace Ipopt
{

/** This class maps a IpoptNLP into one that is used for the restoration phase of %Ipopt.
 *
 *  Given a IpoptNLP
 *  \f{eqnarray*}
 *     \mathrm{min}  && f(x), \\
 *     \mathrm{s.t.} && c(x) = 0,               &\qquad y_c\\
 *                   && d_L \leq d(x) \leq d_U, &\qquad y_d \\
 *                   && x_L \leq  x \leq x_U,   &\qquad z_L, z_U
 *  \f}
 *  parameters \f$\rho\f$, \f$\eta\f$, and a reference point \f$x_r\f$,
 *  a RestoIpoptNLP is the %NLP
 *  \f{eqnarray*}
 *     \mathrm{min}  && \rho (p_c^Te + n_c^Te + p_d^Te + n_d^Te) + \frac{\eta}{2} \Vert D_r (x-x_r) \Vert_2^2, \\
 *     \mathrm{s.t.} && c(x) - p_c + n_c = 0,               &\qquad y_c\\
 *                   && d_L \leq d(x) - p_d + n_d \leq d_U, &\qquad y_d \\
 *                   && x_L \leq  x \leq x_U,               &\qquad z_L, z_U \\
 *                   && p_c, n_c, p_d, n_d \geq 0.          &\qquad ...
 *  \f}
 *  where \f$D_r = \mathrm{diag}(\frac{1}{\max\{1,|x_{r,i}|\}},\ldots,\frac{1}{\max\{1,|x_{r,n}|\}})\f$.
 *
 *  Parameter \f$\rho\f$ is determined by option `resto_penalty_parameter` (default = 1000).
 *  Parameter \f$\eta = \eta_f \sqrt{\mu}\f$, where \f$\eta_f\f$ is determined by option `resto_proximity_weight` (default = 1).
 */
class RestoIpoptNLP: public IpoptNLP
{
public:
   /**@name Constructors/Destructors */
   ///@{
   RestoIpoptNLP(
      IpoptNLP&                  orig_ip_nlp,
      IpoptData&                 orig_ip_data,
      IpoptCalculatedQuantities& orig_ip_cq
   );

   /** Destructor */
   ~RestoIpoptNLP();
   ///@}

   virtual bool Initialize(
      const Journalist&  jnlst,
      const OptionsList& options,
      const std::string& prefix
   );

   /** Initialize (create) structures for the iteration data */
   virtual bool InitializeStructures(
      SmartPtr<Vector>& x,
      bool              init_x,
      SmartPtr<Vector>& y_c,
      bool              init_y_c,
      SmartPtr<Vector>& y_d,
      bool              init_y_d,
      SmartPtr<Vector>& z_L,
      bool              init_z_L,
      SmartPtr<Vector>& z_U,
      bool              init_z_U,
      SmartPtr<Vector>& v_L,
      SmartPtr<Vector>& v_U);

   /** Method accessing the GetWarmStartIterate of the NLP */
   virtual bool GetWarmStartIterate(
      IteratesVector& /*warm_start_iterate*/
   )
   {
      return false;
   }

   void FinalizeSolution(
      SolverReturn               /*status*/,
      const Vector&              /*x*/,
      const Vector&              /*z_L*/,
      const Vector&              /*z_U*/,
      const Vector&              /*c*/,
      const Vector&              /*d*/,
      const Vector&              /*y_c*/,
      const Vector&              /*y_d*/,
      Number                     /*obj_value*/,
      const IpoptData*           /*ip_data*/,
      IpoptCalculatedQuantities* /*ip_cq*/
   )
   { }

   /** Accessor methods for model data */
   ///@{
   /** Method for telling IpoptCalculatedQuantities that the
    *  restoration phase objective function depends on the barrier
    *  parameter
    */
   virtual bool objective_depends_on_mu() const
   {
      return true;
   }

   /** Objective value (incorrect version for restoration phase) */
   virtual Number f(
      const Vector& x
   );

   /** Objective value */
   virtual Number f(
      const Vector& x,
      Number mu
   );

   /** Gradient of the objective (incorrect version for restoration phase) */
   virtual SmartPtr<const Vector> grad_f(
      const Vector& x
   );

   /** Gradient of the objective */
   virtual SmartPtr<const Vector> grad_f(
      const Vector& x,
      Number        mu
   );

   /** Equality constraint residual */
   virtual SmartPtr<const Vector> c(
      const Vector& x
   );

   /** Jacobian Matrix for equality constraints */
   virtual SmartPtr<const Matrix> jac_c(
      const Vector& x
   );

   /** Inequality constraint residual (reformulated
    *  as equalities with slacks */
   virtual SmartPtr<const Vector> d(
      const Vector& x
   );

   /** Jacobian Matrix for inequality constraints */
   virtual SmartPtr<const Matrix> jac_d(
      const Vector& x
   );

   /** Hessian of the Lagrangian (incorrect version for restoration phase) */
   virtual SmartPtr<const SymMatrix> h(
      const Vector& x,
      Number        obj_factor,
      const Vector& yc,
      const Vector& yd
   );

   /** Hessian of the Lagrangian */
   virtual SmartPtr<const SymMatrix> h(
      const Vector& x,
      Number        obj_factor,
      const Vector& yc,
      const Vector& yd,
      Number        mu
   );

   /** Provides a Hessian matrix from the correct matrix space with
    *  uninitialized values.
    *
    *  This can be used in LeastSquareMults to obtain a "zero Hessian".
    */
   virtual SmartPtr<const SymMatrix> uninitialized_h();

   /** Lower bounds on x */
   virtual SmartPtr<const Vector> x_L() const
   {
      return GetRawPtr(x_L_);
   }

   /** Permutation matrix (x_L_ -> x) */
   virtual SmartPtr<const Matrix> Px_L() const
   {
      return GetRawPtr(Px_L_);
   }

   /** Upper bounds on x */
   virtual SmartPtr<const Vector> x_U() const
   {
      return GetRawPtr(x_U_);
   }

   /** Permutation matrix (x_U_ -> x */
   virtual SmartPtr<const Matrix> Px_U() const
   {
      return GetRawPtr(Px_U_);
   }

   /** Lower bounds on d */
   virtual SmartPtr<const Vector> d_L() const
   {
      return GetRawPtr(d_L_);
   }

   /** Permutation matrix (d_L_ -> d) */
   virtual SmartPtr<const Matrix> Pd_L() const
   {
      return GetRawPtr(Pd_L_);
   }

   /** Upper bounds on d */
   virtual SmartPtr<const Vector> d_U() const
   {
      return GetRawPtr(d_U_);
   }

   /** Permutation matrix (d_U_ -> d */
   virtual SmartPtr<const Matrix> Pd_U() const
   {
      return GetRawPtr(Pd_U_);
   }

   virtual SmartPtr<const SymMatrixSpace> HessianMatrixSpace() const
   {
      return GetRawPtr(h_space_);
   }

   virtual SmartPtr<const VectorSpace> x_space() const
   {
      return GetRawPtr(x_space_);
   }
   ///@}

   /** Accessor method for vector/matrix spaces pointers */
   virtual void GetSpaces(
      SmartPtr<const VectorSpace>&    x_space,
      SmartPtr<const VectorSpace>&    c_space,
      SmartPtr<const VectorSpace>&    d_space,
      SmartPtr<const VectorSpace>&    x_l_space,
      SmartPtr<const MatrixSpace>&    px_l_space,
      SmartPtr<const VectorSpace>&    x_u_space,
      SmartPtr<const MatrixSpace>&    px_u_space,
      SmartPtr<const VectorSpace>&    d_l_space,
      SmartPtr<const MatrixSpace>&    pd_l_space,
      SmartPtr<const VectorSpace>&    d_u_space,
      SmartPtr<const MatrixSpace>&    pd_u_space,
      SmartPtr<const MatrixSpace>&    Jac_c_space,
      SmartPtr<const MatrixSpace>&    Jac_d_space,
      SmartPtr<const SymMatrixSpace>& Hess_lagrangian_space);

   /** Method for adapting the variable bounds.
    *
    *  This is called if slacks are becoming too small.
    */
   virtual void AdjustVariableBounds(
      const Vector& new_x_L,
      const Vector& new_x_U,
      const Vector& new_d_L,
      const Vector& new_d_U
   );

   /** User callback method */
   bool IntermediateCallBack(
      AlgorithmMode                       mode,
      Index                               iter,
      Number                              obj_value,
      Number                              inf_pr,
      Number                              inf_du,
      Number                              mu,
      Number                              d_norm,
      Number                              regularization_size,
      Number                              alpha_du,
      Number                              alpha_pr,
      Index                               ls_trials,
      SmartPtr<const IpoptData>           ip_data,
      SmartPtr<IpoptCalculatedQuantities> ip_cq
   );

   /** @name Accessor method for the information of the original NLP.
    *
    *  These methods are not overloaded from IpoptNLP. */
   ///@{
   IpoptNLP& OrigIpNLP() const
   {
      return *orig_ip_nlp_;
   }

   IpoptData& OrigIpData() const
   {
      return *orig_ip_data_;
   }

   IpoptCalculatedQuantities& OrigIpCq() const
   {
      return *orig_ip_cq_;
   }
   ///@}

   /** Accessor Method for obtaining the Rho penalization factor for
    *  the ell_1 norm.
    */
   Number Rho() const
   {
      return rho_;
   }

   /** @name Counters for the number of function evaluations. */
   ///@{
   virtual Index f_evals() const
   {
      return f_evals_;
   }
   virtual Index grad_f_evals() const
   {
      return grad_f_evals_;
   }
   virtual Index c_evals() const
   {
      return c_evals_;
   }
   virtual Index jac_c_evals() const
   {
      return jac_c_evals_;
   }
   virtual Index d_evals() const
   {
      return d_evals_;
   }
   virtual Index jac_d_evals() const
   {
      return jac_d_evals_;
   }
   virtual Index h_evals() const
   {
      return h_evals_;
   }
   ///@}

   /** Method to calculate eta, the factor for the regularization term */
   Number Eta(
      Number mu
   ) const;

   /** Method returning the scaling factors for the 2-norm
    *  penalization term.
    */
   SmartPtr<const Vector> DR_x() const
   {
      return ConstPtr(dr_x_);
   }

   static void RegisterOptions(
      SmartPtr<RegisteredOptions> roptions
   );

private:
   /** @name Pointers for the original NLP information. */
   ///@{
   /** Pointer to the original IpoptNLP */
   SmartPtr<IpoptNLP> orig_ip_nlp_;

   /** Pointer to the original IpoptData */
   SmartPtr<IpoptData> orig_ip_data_;

   /** Pointer to the original IpoptCalculatedQuantities */
   SmartPtr<IpoptCalculatedQuantities> orig_ip_cq_;
   ///@}

   /** Necessary Vector/Matrix spaces */
   ///@{
   SmartPtr<CompoundVectorSpace> x_space_;

   SmartPtr<CompoundVectorSpace> c_space_;

   SmartPtr<CompoundVectorSpace> d_space_;

   SmartPtr<CompoundVectorSpace> x_l_space_;

   SmartPtr<CompoundMatrixSpace> px_l_space_;

   SmartPtr<CompoundVectorSpace> x_u_space_;

   SmartPtr<CompoundMatrixSpace> px_u_space_;

   SmartPtr<CompoundVectorSpace> d_l_space_;

   SmartPtr<CompoundMatrixSpace> pd_l_space_;

   SmartPtr<CompoundVectorSpace> d_u_space_;

   SmartPtr<CompoundMatrixSpace> pd_u_space_;

   SmartPtr<CompoundMatrixSpace> jac_c_space_;

   SmartPtr<CompoundMatrixSpace> jac_d_space_;

   SmartPtr<CompoundSymMatrixSpace> h_space_;
   ///@}

   /**@name Storage for Model Quantities */
   ///@{
   /** Lower bounds on x */
   SmartPtr<CompoundVector> x_L_;

   /** Permutation matrix (x_L_ -> x) */
   SmartPtr<CompoundMatrix> Px_L_;

   /** Upper bounds on x */
   SmartPtr<CompoundVector> x_U_;

   /** Permutation matrix (x_U_ -> x) */
   SmartPtr<CompoundMatrix> Px_U_;

   /** Lower bounds on d */
   SmartPtr<CompoundVector> d_L_;

   /** Permutation matrix (d_L_ -> d) */
   SmartPtr<CompoundMatrix> Pd_L_;

   /** Upper bounds on d */
   SmartPtr<CompoundVector> d_U_;

   /** Permutation matrix (d_U_ -> d */
   SmartPtr<CompoundMatrix> Pd_U_;
   ///@}

   /** @name Values particular to the restoration phase problem statement */
   ///@{
   /** Penalty parameter for the \$l_1\$ norm, given by resto_penalty_parameter */
   Number rho_;

   /** scaling factor for eta calculation */
   Number eta_factor_;

   /** exponent for mu in eta calculation */
   Number eta_mu_exponent_;

   /// Scaling factors for the \$x\$ part of the regularization term
   SmartPtr<const Vector> dr_x_;
   /// Squared scaling factors for the \$x\$ part of the regularization term, for grad_f
   SmartPtr<const Vector> dr2_x_;
   /// Matrix with squared scaling factors, for h()
   SmartPtr<const DiagMatrix> DR2_x_;

   /** \$x\$ part of the reference point in the regularization term */
   SmartPtr<const Vector> x_ref_;
   ///@}

   /**@name Default Compiler Generated Methods
    * (Hidden to avoid implicit creation/calling).
    *
    * These methods are not implemented and
    * we do not want the compiler to implement
    * them for us, so we declare them private
    * and do not define them. This ensures that
    * they will not be implicitly created/called. */
   ///@{
   /** Default Constructor */
   RestoIpoptNLP();

   /** Copy Constructor */
   RestoIpoptNLP(
      const RestoIpoptNLP&
   );

   /** Default Assignment Operator */
   void operator=(
      const RestoIpoptNLP&
   );
   ///@}

   /** @name Algorithmic parameter */
   ///@{
   /** Flag indicating if evaluation of the objective should be
    *  performed for every restoration phase objective function
    *  evaluation.
    */
   bool evaluate_orig_obj_at_resto_trial_;

   /** Flag indicating how Hessian information is obtained */
   HessianApproximationType hessian_approximation_;
   ///@}

   /** Flag indicating if initialization method has been called */
   bool initialized_;

   /** @name Counters for the function evaluations */
   ///@{
   Index f_evals_;
   Index grad_f_evals_;
   Index c_evals_;
   Index jac_c_evals_;
   Index d_evals_;
   Index jac_d_evals_;
   Index h_evals_;
   ///@}
};

} // namespace Ipopt

#endif