ROL
ROL_QuadraticPenalty.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_QUADRATICPENALTY_H
11 #define ROL_QUADRATICPENALTY_H
12 
13 #include "ROL_Objective.hpp"
14 #include "ROL_Constraint.hpp"
15 #include "ROL_Vector.hpp"
16 #include "ROL_Types.hpp"
17 #include "ROL_Ptr.hpp"
18 #include <iostream>
19 
47 namespace ROL {
48 
49 template <class Real>
50 class QuadraticPenalty : public Objective<Real> {
51 private:
52  // Required for quadratic penalty definition
53  const ROL::Ptr<Constraint<Real> > con_;
54  ROL::Ptr<Vector<Real> > multiplier_;
56 
57  // Auxiliary storage
58  ROL::Ptr<Vector<Real> > primalMultiplierVector_;
59  ROL::Ptr<Vector<Real> > dualOptVector_;
60  ROL::Ptr<Vector<Real> > primalConVector_;
61 
62  // Constraint evaluations
63  ROL::Ptr<Vector<Real> > conValue_;
64  Real cscale_;
65 
66  // Evaluation counters
67  int ncval_;
68 
69  // User defined options
70  const bool useScaling_;
71  const int HessianApprox_;
72 
73  // Flags to recompute quantities
75 
76  void evaluateConstraint(const Vector<Real> &x, Real &tol) {
77  if ( !isConstraintComputed_ ) {
78  // Evaluate constraint
79  con_->value(*conValue_,x,tol); ncval_++;
80  isConstraintComputed_ = true;
81  }
82  }
83 
84 public:
85  QuadraticPenalty(const ROL::Ptr<Constraint<Real> > &con,
86  const Vector<Real> &multiplier,
87  const Real penaltyParameter,
88  const Vector<Real> &optVec,
89  const Vector<Real> &conVec,
90  const bool useScaling = false,
91  const int HessianApprox = 0 )
92  : con_(con), penaltyParameter_(penaltyParameter), cscale_(1), ncval_(0),
93  useScaling_(useScaling), HessianApprox_(HessianApprox), isConstraintComputed_(false) {
94 
95  dualOptVector_ = optVec.dual().clone();
96  primalConVector_ = conVec.clone();
97  conValue_ = conVec.clone();
98  multiplier_ = multiplier.clone();
99  primalMultiplierVector_ = multiplier.clone();
100  }
101 
102  void setScaling(const Real cscale = 1) {
103  cscale_ = cscale;
104  }
105 
106  virtual void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
107  con_->update(x,flag,iter);
108  isConstraintComputed_ = ((flag || (!flag && iter < 0)) ? false : isConstraintComputed_);
109  }
110 
111  virtual Real value( const Vector<Real> &x, Real &tol ) {
112  // Evaluate constraint
113  evaluateConstraint(x,tol);
114  // Apply Lagrange multiplier to constraint
115  Real cval = cscale_*multiplier_->dot(conValue_->dual());
116  // Compute penalty term
117  Real pval = cscale_*cscale_*conValue_->dot(*conValue_);
118  // Compute quadratic penalty value
119  const Real half(0.5);
120  Real val(0);
121  if (useScaling_) {
122  val = cval/penaltyParameter_ + half*pval;
123  }
124  else {
125  val = cval + half*penaltyParameter_*pval;
126  }
127  return val;
128  }
129 
130  virtual void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
131  // Evaluate constraint
132  evaluateConstraint(x,tol);
133  // Compute gradient of Augmented Lagrangian
134  primalMultiplierVector_->set(conValue_->dual());
135  if ( useScaling_ ) {
138  }
139  else {
142  }
143  con_->applyAdjointJacobian(g,*primalMultiplierVector_,x,tol);
144  }
145 
146  virtual void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
147  // Apply objective Hessian to a vector
148  if (HessianApprox_ < 3) {
149  con_->update(x);
150  con_->applyJacobian(*primalConVector_,v,x,tol);
151  con_->applyAdjointJacobian(hv,primalConVector_->dual(),x,tol);
152  if (!useScaling_) {
154  }
155  else {
156  hv.scale(cscale_*cscale_);
157  }
158 
159  if (HessianApprox_ == 1) {
160  // Apply Augmented Lagrangian Hessian to a vector
162  if ( useScaling_ ) {
164  }
165  else {
167  }
168  con_->applyAdjointHessian(*dualOptVector_,*primalMultiplierVector_,v,x,tol);
169  hv.plus(*dualOptVector_);
170  }
171 
172  if (HessianApprox_ == 0) {
173  // Evaluate constraint
174  evaluateConstraint(x,tol);
175  // Apply Augmented Lagrangian Hessian to a vector
176  primalMultiplierVector_->set(conValue_->dual());
177  if ( useScaling_ ) {
180  }
181  else {
184  }
185  con_->applyAdjointHessian(*dualOptVector_,*primalMultiplierVector_,v,x,tol);
186  hv.plus(*dualOptVector_);
187  }
188  }
189  else {
190  hv.zero();
191  }
192  }
193 
194  // Return constraint value
195  virtual void getConstraintVec(Vector<Real> &c, const Vector<Real> &x) {
196  Real tol = std::sqrt(ROL_EPSILON<Real>());
197  // Evaluate constraint
198  evaluateConstraint(x,tol);
199  c.set(*conValue_);
200  }
201 
202  // Return total number of constraint evaluations
203  virtual int getNumberConstraintEvaluations(void) const {
204  return ncval_;
205  }
206 
207  // Reset with upated penalty parameter
208  virtual void reset(const Vector<Real> &multiplier, const Real penaltyParameter) {
209  ncval_ = 0;
210  multiplier_->set(multiplier);
211  penaltyParameter_ = penaltyParameter;
212  }
213 }; // class AugmentedLagrangian
214 
215 } // namespace ROL
216 
217 #endif
Provides the interface to evaluate objective functions.
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:192
virtual void scale(const Real alpha)=0
Compute where .
virtual void reset(const Vector< Real > &multiplier, const Real penaltyParameter)
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
virtual void plus(const Vector &x)=0
Compute , where .
ROL::Ptr< Vector< Real > > primalConVector_
QuadraticPenalty(const ROL::Ptr< Constraint< Real > > &con, const Vector< Real > &multiplier, const Real penaltyParameter, const Vector< Real > &optVec, const Vector< Real > &conVec, const bool useScaling=false, const int HessianApprox=0)
Provides the interface to evaluate the quadratic constraint penalty.
Contains definitions of custom data types in ROL.
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
const ROL::Ptr< Constraint< Real > > con_
void evaluateConstraint(const Vector< Real > &x, Real &tol)
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
ROL::Ptr< Vector< Real > > multiplier_
virtual int getNumberConstraintEvaluations(void) const
virtual Real value(const Vector< Real > &x, Real &tol)
Compute value.
void setScaling(const Real cscale=1)
ROL::Ptr< Vector< Real > > conValue_
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
virtual void getConstraintVec(Vector< Real > &c, const Vector< Real > &x)
ROL::Ptr< Vector< Real > > primalMultiplierVector_
Defines the general constraint operator interface.
ROL::Ptr< Vector< Real > > dualOptVector_