ROL
ROL_DynamicConstraint.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 #pragma once
11 #ifndef ROL_DYNAMICCONSTRAINT_HPP
12 #define ROL_DYNAMICCONSTRAINT_HPP
13 
14 
16 #include "ROL_Types.hpp"
17 #include "ROL_DynamicFunction.hpp"
20 
21 namespace ROL {
22 
23 template<typename Real>
25 
26 template<typename Real>
28 }
29 
32 
50 namespace ROL {
51 
52 template<typename Real>
53 class DynamicConstraint : public DynamicFunction<Real> {
54 private:
55  // Additional vector storage for solve
56  Ptr<Vector<Real>> unew_;
57  Ptr<Vector<Real>> jv_;
58 
59  // Default parameters for solve (backtracking Newton)
60  const Real DEFAULT_atol_;
61  const Real DEFAULT_rtol_;
62  const Real DEFAULT_stol_;
63  const Real DEFAULT_factor_;
64  const Real DEFAULT_decr_;
65  const int DEFAULT_maxit_;
66  const bool DEFAULT_print_;
67  const bool DEFAULT_zero_;
69 
70  // User-set parameters for solve (backtracking Newton)
71  Real atol_;
72  Real rtol_;
73  Real stol_;
74  Real factor_;
75  Real decr_;
76  int maxit_;
77  bool print_;
78  bool zero_;
80 
81  // Flag to initialize vector storage in solve
83 
84 public:
85 
86  using V = Vector<Real>;
89 
90  virtual ~DynamicConstraint() {}
91 
92  DynamicConstraint( std::initializer_list<std::string> zero_deriv_terms={} ):
93  DynamicFunction<Real>(zero_deriv_terms),
94  unew_ ( nullPtr ),
95  jv_ ( nullPtr ),
96  DEFAULT_atol_ ( 1e-4*std::sqrt(ROL_EPSILON<Real>()) ),
97  DEFAULT_rtol_ ( 1e0 ),
98  DEFAULT_stol_ ( std::sqrt(ROL_EPSILON<Real>()) ),
99  DEFAULT_factor_ ( 0.5 ),
100  DEFAULT_decr_ ( 1e-4 ),
101  DEFAULT_maxit_ ( 500 ),
102  DEFAULT_print_ ( false ),
103  DEFAULT_zero_ ( false ),
104  DEFAULT_solverType_ ( 0 ),
105  atol_ ( DEFAULT_atol_ ),
106  rtol_ ( DEFAULT_rtol_ ),
107  stol_ ( DEFAULT_stol_ ),
109  decr_ ( DEFAULT_decr_ ),
112  zero_ ( DEFAULT_zero_ ),
114  firstSolve_ ( true ) {}
115 
116 
117 
118  virtual void update( const V& uo, const V& un, const V& z, const TS& ts ) {
119  update_uo( uo, ts );
120  update_un( un, ts );
121  update_z( z, ts );
122  }
123 
127 
128  virtual void value( V& c, const V& uo, const V& un,
129  const V& z, const TS& ts ) const = 0;
130 
131  virtual void solve( V& c, const V& uo, V& un,
132  const V& z, const TS& ts ) {
133  if ( zero_ ) un.zero();
134  Ptr<std::ostream> stream = makeStreamPtr(std::cout, print_);
135  update(uo,un,z,ts);
136  value(c,uo,un,z,ts);
137  Real cnorm = c.norm();
138  const Real ctol = std::min(atol_, rtol_*cnorm);
139  if (solverType_==0 || solverType_==3 || solverType_==4) {
140  if ( firstSolve_ ) {
141  unew_ = un.clone();
142  jv_ = un.clone();
143  firstSolve_ = false;
144  }
145  Real alpha(1), tmp(0);
146  int cnt = 0;
147  if ( print_ ) {
148  *stream << " Default DynamicConstraint::solve" << std::endl;
149  *stream << " ";
150  *stream << std::setw(6) << std::left << "iter";
151  *stream << std::setw(15) << std::left << "rnorm";
152  *stream << std::setw(15) << std::left << "alpha";
153  *stream << std::endl;
154  }
155  for (cnt = 0; cnt < maxit_; ++cnt) {
156  // Compute Newton step
157  applyInverseJacobian_un(*jv_,c,uo,un,z,ts);
158  unew_->set(un);
159  unew_->axpy(-alpha, *jv_);
160  //update_un(*unew_,ts);
161  update(uo,*unew_,z,ts);
162  value(c,uo,*unew_,z,ts);
163  tmp = c.norm();
164  // Perform backtracking line search
165  while ( tmp > (1.0-decr_*alpha)*cnorm &&
166  alpha > stol_ ) {
167  alpha *= factor_;
168  unew_->set(un);
169  unew_->axpy(-alpha,*jv_);
170  //update_un(*unew_,ts);
171  update(uo,*unew_,z,ts);
172  value(c,uo,*unew_,z,ts);
173  tmp = c.norm();
174  }
175  if ( print_ ) {
176  *stream << " ";
177  *stream << std::setw(6) << std::left << cnt;
178  *stream << std::scientific << std::setprecision(6);
179  *stream << std::setw(15) << std::left << tmp;
180  *stream << std::scientific << std::setprecision(6);
181  *stream << std::setw(15) << std::left << alpha;
182  *stream << std::endl;
183  }
184  // Update iterate
185  cnorm = tmp;
186  un.set(*unew_);
187  if (cnorm <= ctol) { // = covers the case of identically zero residual
188  break;
189  }
190  update(uo,un,z,ts);
191  alpha = 1.0;
192  }
193  }
194  if (solverType_==1 || (solverType_==3 && cnorm > ctol)) {
195  Ptr<DynamicConstraint<Real>> con = makePtrFromRef(*this);
196  Ptr<Objective<Real>> obj
197  = makePtr<NonlinearLeastSquaresObjective_Dynamic<Real>>(con,c,makePtrFromRef(uo),makePtrFromRef(z),makePtrFromRef(ts),true);
198  ParameterList parlist;
199  parlist.sublist("General").set("Output Level",1);
200  parlist.sublist("General").sublist("Krylov").set("Iteration Limit",100);
201  parlist.sublist("Step").sublist("Trust Region").set("Subproblem Solver","Truncated CG");
202  parlist.sublist("Status Test").set("Gradient Tolerance",ctol);
203  parlist.sublist("Status Test").set("Step Tolerance",stol_);
204  parlist.sublist("Status Test").set("Iteration Limit",maxit_);
205  Ptr<TypeU::Algorithm<Real>> algo = makePtr<TypeU::TrustRegionAlgorithm<Real>>(parlist);
206  algo->run(un,*obj,*stream);
207  value(c,uo,un,z,ts);
208  }
209  if (solverType_==2 || (solverType_==4 && cnorm > ctol)) {
210  Ptr<Constraint<Real>> con
211  = makePtr<Constraint_DynamicState<Real>>(makePtrFromRef(*this),makePtrFromRef(uo),makePtrFromRef(z),makePtrFromRef(ts));
212  Ptr<Objective<Real>> obj = makePtr<Objective_FSsolver<Real>>();
213  Ptr<Vector<Real>> l = c.dual().clone();
214  ParameterList parlist;
215  parlist.sublist("General").set("Output Level",1);
216  parlist.sublist("Status Test").set("Constraint Tolerance",ctol);
217  parlist.sublist("Status Test").set("Step Tolerance",stol_);
218  parlist.sublist("Status Test").set("Iteration Limit",maxit_);
219  Ptr<TypeE::Algorithm<Real>> algo = makePtr<TypeE::AugmentedLagrangianAlgorithm<Real>>(parlist);
220  algo->run(un,*obj,*con,*l,*stream);
221  value(c,uo,un,z,ts);
222  }
223  if (solverType_ > 4 || solverType_ < 0) {
224  ROL_TEST_FOR_EXCEPTION(true, std::invalid_argument,
225  ">>> ERROR (ROL:DynamicConstraint:solve): Invalid solver type!");
226  }
227  }
228 
249  virtual void setSolveParameters(ParameterList &parlist) {
250  ParameterList & list = parlist.sublist("Dynamic Constraint").sublist("Solve");
251  atol_ = list.get("Absolute Residual Tolerance", DEFAULT_atol_);
252  rtol_ = list.get("Relative Residual Tolerance", DEFAULT_rtol_);
253  maxit_ = list.get("Iteration Limit", DEFAULT_maxit_);
254  decr_ = list.get("Sufficient Decrease Tolerance", DEFAULT_decr_);
255  stol_ = list.get("Step Tolerance", DEFAULT_stol_);
256  factor_ = list.get("Backtracking Factor", DEFAULT_factor_);
257  print_ = list.get("Output Iteration History", DEFAULT_print_);
258  zero_ = list.get("Zero Initial Guess", DEFAULT_zero_);
259  solverType_ = list.get("Solver Type", DEFAULT_solverType_);
260  }
261 
262  //----------------------------------------------------------------------------
263  // Partial Jacobians
264  virtual void applyJacobian_uo( V& jv, const V& vo, const V& uo,
265  const V& un, const V& z,
266  const TS& ts ) const {}
267 
268  virtual void applyJacobian_un( V& jv, const V& vn, const V& uo,
269  const V& un, const V& z,
270  const TS& ts ) const {}
271 
272  virtual void applyJacobian_z( V& jv, const V& vz, const V& uo,
273  const V& un, const V& z,
274  const TS& ts ) const {}
275 
276  //----------------------------------------------------------------------------
277  // Adjoint partial Jacobians
278  virtual void applyAdjointJacobian_uo( V& ajv, const V& dualv, const V& uo,
279  const V& un, const V& z,
280  const TS& ts ) const {}
281 
282  virtual void applyAdjointJacobian_un( V& ajv, const V& dualv, const V& uo,
283  const V& un, const V& z,
284  const TS& ts ) const {}
285 
286  virtual void applyAdjointJacobian_z( V& ajv, const V& dualv, const V& uo,
287  const V& un, const V& z,
288  const TS& ts ) const {}
289 
290  //----------------------------------------------------------------------------
291  // Inverses
292  virtual void applyInverseJacobian_un( V& ijv, const V& vn, const V& uo,
293  const V& un, const V& z,
294  const TS& ts ) const {}
295 
296  virtual void applyInverseAdjointJacobian_un( V& iajv, const V& vn, const V& uo,
297  const V& un, const V& z,
298  const TS& ts ) const {}
299 
300  //----------------------------------------------------------------------------
301  // Adjoint Hessian components
302  virtual void applyAdjointHessian_un_un( V& ahwv, const V& wn, const V& vn,
303  const V& uo, const V& un,
304  const V& z, const TS& ts ) const {
305  ahwv.zero();
306  }
307 
308  virtual void applyAdjointHessian_un_uo( V& ahwv, const V& w, const V& vn,
309  const V& uo, const V& un,
310  const V& z, const TS& ts ) const {
311  ahwv.zero();
312  }
313 
314  virtual void applyAdjointHessian_un_z( V& ahwv, const V& w, const V& vn,
315  const V& uo, const V& un,
316  const V& z, const TS& ts ) const {
317  ahwv.zero();
318  }
319 
320  virtual void applyAdjointHessian_uo_un( V& ahwv, const V& w, const V& vo,
321  const V& uo, const V& un,
322  const V& z, const TS& ts ) const {
323  ahwv.zero();
324  }
325 
326  virtual void applyAdjointHessian_uo_uo( V& ahwv, const V& w, const V& v,
327  const V& uo, const V& un,
328  const V& z, const TS& ts ) const {
329  ahwv.zero();
330  }
331 
332  virtual void applyAdjointHessian_uo_z( V& ahwv, const V& w, const V& vo,
333  const V& uo, const V& un,
334  const V& z, const TS& ts ) const {
335  ahwv.zero();
336  }
337 
338  virtual void applyAdjointHessian_z_un( V& ahwv, const V& w, const V& vz,
339  const V& uo, const V& un,
340  const V& z, const TS& ts ) const {
341  ahwv.zero();
342  }
343 
344  virtual void applyAdjointHessian_z_uo( V& ahwv, const V& w, const V& vz,
345  const V& uo, const V& un,
346  const V& z, const TS& ts ) const {
347  ahwv.zero();
348  }
349 
350  virtual void applyAdjointHessian_z_z( V& ahwv, const V& w, const V& vz,
351  const V& uo, const V& un,
352  const V& z, const TS& ts ) const {
353  ahwv.zero();
354  }
355 
356 }; // DynamicConstraint
357 
358 
359 } // namespace ROL
360 
361 
362 #endif // ROL_DYNAMICCONSTRAINT_HPP
363 
Ptr< Vector< Real > > jv_
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 applyAdjointHessian_un_z(V &ahwv, const V &w, const V &vn, const V &uo, const V &un, const V &z, const TS &ts) const
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void solve(V &c, const V &uo, V &un, const V &z, const TS &ts)
Defines the time-dependent constraint operator interface for simulation-based optimization.
virtual void applyAdjointHessian_uo_un(V &ahwv, const V &w, const V &vo, const V &uo, const V &un, const V &z, const TS &ts) const
Defines the linear algebra of vector space on a generic partitioned vector.
virtual void applyAdjointHessian_z_z(V &ahwv, const V &w, const V &vz, const V &uo, const V &un, const V &z, const TS &ts) const
Provides update interface, casting and vector management to DynamicConstraint and DynamicObjective...
Ptr< Vector< Real > > unew_
Contains definitions of custom data types in ROL.
virtual void applyAdjointHessian_z_un(V &ahwv, const V &w, const V &vz, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void applyAdjointJacobian_un(V &ajv, const V &dualv, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void applyJacobian_uo(V &jv, const V &vo, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void update(const V &uo, const V &un, const V &z, const TS &ts)
virtual void applyInverseAdjointJacobian_un(V &iajv, const V &vn, const V &uo, const V &un, const V &z, const TS &ts) const
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
Ptr< ostream > makeStreamPtr(ostream &os, bool noSuppressOutput=true)
Definition: ROL_Stream.hpp:41
virtual void applyAdjointJacobian_z(V &ajv, const V &dualv, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void setSolveParameters(ParameterList &parlist)
Set solve parameters.
virtual void applyInverseJacobian_un(V &ijv, const V &vn, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void update_z(const V &x, const TS &ts)
virtual void value(V &c, const V &uo, const V &un, const V &z, const TS &ts) const =0
virtual void update_uo(const V &x, const TS &ts)
DynamicConstraint(std::initializer_list< std::string > zero_deriv_terms={})
virtual void applyAdjointJacobian_uo(V &ajv, const V &dualv, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void applyJacobian_un(V &jv, const V &vn, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void applyJacobian_z(V &jv, const V &vz, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void applyAdjointHessian_un_uo(V &ahwv, const V &w, const V &vn, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void applyAdjointHessian_z_uo(V &ahwv, const V &w, const V &vz, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void applyAdjointHessian_uo_uo(V &ahwv, const V &w, const V &v, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
virtual Real norm() const =0
Returns where .
virtual void applyAdjointHessian_un_un(V &ahwv, const V &wn, const V &vn, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void applyAdjointHessian_uo_z(V &ahwv, const V &w, const V &vo, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void update_un(const V &x, const TS &ts)