ROL
ROL_QuadraticPenalty_SimOpt.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_SIMOPT_H
11 #define ROL_QUADRATICPENALTY_SIMOPT_H
12 
13 #include "ROL_Objective_SimOpt.hpp"
15 #include "ROL_Vector.hpp"
16 #include "ROL_Types.hpp"
17 #include "ROL_Ptr.hpp"
18 #include <iostream>
19 
54 namespace ROL {
55 
56 template <class Real>
58 private:
59  // Required for quadratic penalty definition
60  const ROL::Ptr<Constraint_SimOpt<Real> > con_;
61  ROL::Ptr<Vector<Real> > multiplier_;
63 
64  // Auxiliary storage
65  ROL::Ptr<Vector<Real> > primalMultiplierVector_;
66  ROL::Ptr<Vector<Real> > dualSimVector_;
67  ROL::Ptr<Vector<Real> > dualOptVector_;
68  ROL::Ptr<Vector<Real> > primalConVector_;
69 
70  // Constraint evaluations
71  ROL::Ptr<Vector<Real> > conValue_;
72 
73  // Evaluation counters
74  int ncval_;
75 
76  // User defined options
77  const bool useScaling_;
78  const int HessianApprox_;
79 
80  // Flags to recompute quantities
82 
83  void evaluateConstraint(const Vector<Real> &u, const Vector<Real> &z, Real &tol) {
84  if ( !isConstraintComputed_ ) {
85  // Evaluate constraint
86  con_->value(*conValue_,u,z,tol); ncval_++;
87  isConstraintComputed_ = true;
88  }
89  }
90 
91 public:
93  const Vector<Real> &multiplier,
94  const Real penaltyParameter,
95  const Vector<Real> &simVec,
96  const Vector<Real> &optVec,
97  const Vector<Real> &conVec,
98  const bool useScaling = false,
99  const int HessianApprox = 0 )
100  : con_(con), penaltyParameter_(penaltyParameter), ncval_(0),
101  useScaling_(useScaling), HessianApprox_(HessianApprox), isConstraintComputed_(false) {
102 
103  dualSimVector_ = simVec.dual().clone();
104  dualOptVector_ = optVec.dual().clone();
105  primalConVector_ = conVec.clone();
106  conValue_ = conVec.clone();
107  multiplier_ = multiplier.clone();
108  primalMultiplierVector_ = multiplier.clone();
109  }
110 
111  virtual void update( const Vector<Real> &u, const Vector<Real> &z, bool flag = true, int iter = -1 ) {
112  con_->update(u,z,flag,iter);
113  isConstraintComputed_ = ( flag ? false : isConstraintComputed_ );
114  }
115 
116  virtual Real value( const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
117  // Evaluate constraint
118  evaluateConstraint(u,z,tol);
119  // Apply Lagrange multiplier to constraint
120  Real cval = multiplier_->dot(conValue_->dual());
121  // Compute penalty term
122  Real pval = conValue_->dot(*conValue_);
123  // Compute quadratic penalty value
124  const Real half(0.5);
125  Real val(0);
126  if (useScaling_) {
127  val = cval/penaltyParameter_ + half*pval;
128  }
129  else {
130  val = cval + half*penaltyParameter_*pval;
131  }
132  return val;
133  }
134 
135  virtual void gradient_1( Vector<Real> &g, const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
136  // Evaluate constraint
137  evaluateConstraint(u,z,tol);
138  // Compute gradient of Augmented Lagrangian
139  primalMultiplierVector_->set(conValue_->dual());
140  if ( useScaling_ ) {
141  primalMultiplierVector_->axpy(static_cast<Real>(1)/penaltyParameter_,*multiplier_);
142  }
143  else {
146  }
147  con_->applyAdjointJacobian_1(g,*primalMultiplierVector_,u,z,tol);
148  }
149 
150  virtual void gradient_2( Vector<Real> &g, const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
151  // Evaluate constraint
152  evaluateConstraint(u,z,tol);
153  // Compute gradient of Augmented Lagrangian
154  primalMultiplierVector_->set(conValue_->dual());
155  if ( useScaling_ ) {
156  primalMultiplierVector_->axpy(static_cast<Real>(1)/penaltyParameter_,*multiplier_);
157  }
158  else {
161  }
162  con_->applyAdjointJacobian_2(g,*primalMultiplierVector_,u,z,tol);
163  }
164 
165  virtual void hessVec_11( Vector<Real> &hv, const Vector<Real> &v,
166  const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
167  // Apply objective Hessian to a vector
168  if (HessianApprox_ < 2) {
169  con_->applyJacobian_1(*primalConVector_,v,u,z,tol);
170  con_->applyAdjointJacobian_1(hv,primalConVector_->dual(),u,z,tol);
171  if (!useScaling_) {
173  }
174 
175  if (HessianApprox_ == 0) {
176  // Evaluate constraint
177  evaluateConstraint(u,z,tol);
178  // Apply Augmented Lagrangian Hessian to a vector
179  primalMultiplierVector_->set(conValue_->dual());
180  if ( useScaling_ ) {
181  primalMultiplierVector_->axpy(static_cast<Real>(1)/penaltyParameter_,*multiplier_);
182  }
183  else {
186  }
187  con_->applyAdjointHessian_11(*dualSimVector_,*primalMultiplierVector_,v,u,z,tol);
188  hv.plus(*dualSimVector_);
189  }
190  }
191  else {
192  hv.zero();
193  }
194  }
195 
196  virtual void hessVec_12( Vector<Real> &hv, const Vector<Real> &v,
197  const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
198  // Apply objective Hessian to a vector
199  if (HessianApprox_ < 2) {
200  con_->applyJacobian_2(*primalConVector_,v,u,z,tol);
201  con_->applyAdjointJacobian_1(hv,primalConVector_->dual(),u,z,tol);
202  if (!useScaling_) {
204  }
205 
206  if (HessianApprox_ == 0) {
207  // Evaluate constraint
208  evaluateConstraint(u,z,tol);
209  // Apply Augmented Lagrangian Hessian to a vector
210  primalMultiplierVector_->set(conValue_->dual());
211  if ( useScaling_ ) {
212  primalMultiplierVector_->axpy(static_cast<Real>(1)/penaltyParameter_,*multiplier_);
213  }
214  else {
217  }
218  con_->applyAdjointHessian_21(*dualSimVector_,*primalMultiplierVector_,v,u,z,tol);
219  hv.plus(*dualSimVector_);
220  }
221  }
222  else {
223  hv.zero();
224  }
225  }
226 
227  virtual void hessVec_21( Vector<Real> &hv, const Vector<Real> &v,
228  const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
229  // Apply objective Hessian to a vector
230  if (HessianApprox_ < 2) {
231  con_->applyJacobian_1(*primalConVector_,v,u,z,tol);
232  con_->applyAdjointJacobian_2(hv,primalConVector_->dual(),u,z,tol);
233  if (!useScaling_) {
235  }
236 
237  if (HessianApprox_ == 0) {
238  // Evaluate constraint
239  evaluateConstraint(u,z,tol);
240  // Apply Augmented Lagrangian Hessian to a vector
241  primalMultiplierVector_->set(conValue_->dual());
242  if ( useScaling_ ) {
243  primalMultiplierVector_->axpy(static_cast<Real>(1)/penaltyParameter_,*multiplier_);
244  }
245  else {
248  }
249  con_->applyAdjointHessian_12(*dualOptVector_,*primalMultiplierVector_,v,u,z,tol);
250  hv.plus(*dualOptVector_);
251  }
252  }
253  else {
254  hv.zero();
255  }
256  }
257 
258  virtual void hessVec_22( Vector<Real> &hv, const Vector<Real> &v,
259  const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
260  // Apply objective Hessian to a vector
261  if (HessianApprox_ < 2) {
262  con_->applyJacobian_2(*primalConVector_,v,u,z,tol);
263  con_->applyAdjointJacobian_2(hv,primalConVector_->dual(),u,z,tol);
264  if (!useScaling_) {
266  }
267 
268  if (HessianApprox_ == 0) {
269  // Evaluate constraint
270  evaluateConstraint(u,z,tol);
271  // Apply Augmented Lagrangian Hessian to a vector
272  primalMultiplierVector_->set(conValue_->dual());
273  if ( useScaling_ ) {
274  primalMultiplierVector_->axpy(static_cast<Real>(1)/penaltyParameter_,*multiplier_);
275  }
276  else {
279  }
280  con_->applyAdjointHessian_22(*dualOptVector_,*primalMultiplierVector_,v,u,z,tol);
281  hv.plus(*dualOptVector_);
282  }
283  }
284  else {
285  hv.zero();
286  }
287  }
288 
289  // Return constraint value
290  virtual void getConstraintVec(Vector<Real> &c, const Vector<Real> &u, const Vector<Real> &z) {
291  Real tol = std::sqrt(ROL_EPSILON<Real>());
292  // Evaluate constraint
293  evaluateConstraint(u,z,tol);
294  c.set(*conValue_);
295  }
296 
297  // Return total number of constraint evaluations
298  virtual int getNumberConstraintEvaluations(void) const {
299  return ncval_;
300  }
301 
302  // Reset with upated penalty parameter
303  virtual void reset(const Vector<Real> &multiplier, const Real penaltyParameter) {
304  ncval_ = 0;
305  multiplier_->set(multiplier);
306  penaltyParameter_ = penaltyParameter;
307  }
308 }; // class AugmentedLagrangian
309 
310 } // namespace ROL
311 
312 #endif
Provides the interface to evaluate simulation-based 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
const ROL::Ptr< Constraint_SimOpt< Real > > con_
virtual void scale(const Real alpha)=0
Compute where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void hessVec_21(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
virtual void plus(const Vector &x)=0
Compute , where .
Contains definitions of custom data types in ROL.
virtual void hessVec_22(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Provides the interface to evaluate the quadratic SimOpt constraint penalty.
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
ROL::Ptr< Vector< Real > > multiplier_
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual int getNumberConstraintEvaluations(void) const
ROL::Ptr< Vector< Real > > primalConVector_
ROL::Ptr< Vector< Real > > primalMultiplierVector_
virtual void hessVec_11(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply Hessian approximation to vector.
ROL::Ptr< Vector< Real > > dualOptVector_
virtual Real value(const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Compute value.
void evaluateConstraint(const Vector< Real > &u, const Vector< Real > &z, Real &tol)
virtual void hessVec_12(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
virtual void getConstraintVec(Vector< Real > &c, const Vector< Real > &u, const Vector< Real > &z)
ROL::Ptr< Vector< Real > > dualSimVector_
virtual void reset(const Vector< Real > &multiplier, const Real penaltyParameter)
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1)
Update objective function. u is an iterate, z is an iterate, flag = true if the iterate has changed...
Defines the constraint operator interface for simulation-based optimization.
virtual void gradient_1(Vector< Real > &g, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Compute gradient with respect to first component.
virtual void gradient_2(Vector< Real > &g, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Compute gradient with respect to second component.
QuadraticPenalty_SimOpt(const ROL::Ptr< Constraint_SimOpt< Real > > &con, const Vector< Real > &multiplier, const Real penaltyParameter, const Vector< Real > &simVec, const Vector< Real > &optVec, const Vector< Real > &conVec, const bool useScaling=false, const int HessianApprox=0)