ROL
ROL_DynamicConstraint.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
44 
45 #pragma once
46 #ifndef ROL_DYNAMICCONSTRAINT_HPP
47 #define ROL_DYNAMICCONSTRAINT_HPP
48 
49 #include "ROL_DynamicFunction.hpp"
50 #include "ROL_TimeStamp.hpp"
51 
55 #include "ROL_Algorithm.hpp"
56 #include "ROL_Types.hpp"
57 
76 namespace ROL {
77 
78 template<typename Real>
79 class DynamicConstraint : public DynamicFunction<Real> {
80 private:
81  // Additional vector storage for solve
82  Ptr<Vector<Real>> unew_;
83  Ptr<Vector<Real>> jv_;
84 
85  // Default parameters for solve (backtracking Newton)
86  const Real DEFAULT_atol_;
87  const Real DEFAULT_rtol_;
88  const Real DEFAULT_stol_;
89  const Real DEFAULT_factor_;
90  const Real DEFAULT_decr_;
91  const int DEFAULT_maxit_;
92  const bool DEFAULT_print_;
93  const bool DEFAULT_zero_;
95 
96  // User-set parameters for solve (backtracking Newton)
97  Real atol_;
98  Real rtol_;
99  Real stol_;
100  Real factor_;
101  Real decr_;
102  int maxit_;
103  bool print_;
104  bool zero_;
106 
107  // Flag to initialize vector storage in solve
109 
110 public:
111 
112  using V = Vector<Real>;
115 
116  virtual ~DynamicConstraint() {}
117 
119  : unew_ ( nullPtr ),
120  jv_ ( nullPtr ),
121  DEFAULT_atol_ ( 1e-4*std::sqrt(ROL_EPSILON<Real>()) ),
122  DEFAULT_rtol_ ( 1e0 ),
123  DEFAULT_stol_ ( std::sqrt(ROL_EPSILON<Real>()) ),
124  DEFAULT_factor_ ( 0.5 ),
125  DEFAULT_decr_ ( 1e-4 ),
126  DEFAULT_maxit_ ( 500 ),
127  DEFAULT_print_ ( false ),
128  DEFAULT_zero_ ( false ),
129  DEFAULT_solverType_ ( 0 ),
130  atol_ ( DEFAULT_atol_ ),
131  rtol_ ( DEFAULT_rtol_ ),
132  stol_ ( DEFAULT_stol_ ),
134  decr_ ( DEFAULT_decr_ ),
137  zero_ ( DEFAULT_zero_ ),
139  firstSolve_ ( true ) {}
140 
141  virtual void update( const V& uo, const V& un, const V& z, const TS& ts ) {
142  update_uo( uo, ts );
143  update_un( un, ts );
144  update_z( z, ts );
145  }
146 
147  virtual void update_uo( const V& uo, const TS& ts ) { }
148  virtual void update_un( const V& un, const TS& ts ) { }
149  virtual void update_z( const V& z, const TS& ts ) { }
150 
151  virtual void value( V& c, const V& uo, const V& un,
152  const V& z, const TS& ts ) const = 0;
153 
154  virtual void solve( V& c, const V& uo, V& un,
155  const V& z, const TS& ts ) {
156  if ( zero_ ) {
157  un.zero();
158  }
159  update(uo,un,z,ts);
160  value(c,uo,un,z,ts);
161  Real cnorm = c.norm();
162  const Real ctol = std::min(atol_, rtol_*cnorm);
163  if (solverType_==0 || solverType_==3 || solverType_==4) {
164  if ( firstSolve_ ) {
165  unew_ = un.clone();
166  jv_ = un.clone();
167  firstSolve_ = false;
168  }
169  Real alpha(1), tmp(0);
170  int cnt = 0;
171  if ( print_ ) {
172  std::cout << "\n Default DynamicConstraint::solve\n";
173  std::cout << " ";
174  std::cout << std::setw(6) << std::left << "iter";
175  std::cout << std::setw(15) << std::left << "rnorm";
176  std::cout << std::setw(15) << std::left << "alpha";
177  std::cout << "\n";
178  }
179  for (cnt = 0; cnt < maxit_; ++cnt) {
180  // Compute Newton step
181  applyInverseJacobian_un(*jv_,c,uo,un,z,ts);
182  unew_->set(un);
183  unew_->axpy(-alpha, *jv_);
184  //update_un(*unew_,ts);
185  update(uo,*unew_,z,ts);
186  value(c,uo,*unew_,z,ts);
187  tmp = c.norm();
188  // Perform backtracking line search
189  while ( tmp > (1.0-decr_*alpha)*cnorm &&
190  alpha > stol_ ) {
191  alpha *= factor_;
192  unew_->set(un);
193  unew_->axpy(-alpha,*jv_);
194  //update_un(*unew_,ts);
195  update(uo,*unew_,z,ts);
196  value(c,uo,*unew_,z,ts);
197  tmp = c.norm();
198  }
199  if ( print_ ) {
200  std::cout << " ";
201  std::cout << std::setw(6) << std::left << cnt;
202  std::cout << std::scientific << std::setprecision(6);
203  std::cout << std::setw(15) << std::left << tmp;
204  std::cout << std::scientific << std::setprecision(6);
205  std::cout << std::setw(15) << std::left << alpha;
206  std::cout << "\n";
207  }
208  // Update iterate
209  cnorm = tmp;
210  un.set(*unew_);
211  if (cnorm <= ctol) { // = covers the case of identically zero residual
212  break;
213  }
214  update(uo,un,z,ts);
215  alpha = 1.0;
216  }
217  }
218  if (solverType_==1 || (solverType_==3 && cnorm > ctol)) {
219  ROL::Ptr<DynamicConstraint<Real>> con_ptr = ROL::makePtrFromRef(*this);
220  ROL::Ptr<const Vector<Real>> uo_ptr = ROL::makePtrFromRef(uo);
221  ROL::Ptr<const Vector<Real>> z_ptr = ROL::makePtrFromRef(z);
222  ROL::Ptr<const TimeStamp<Real>> ts_ptr = ROL::makePtrFromRef(ts);
223  ROL::Ptr<Objective<Real>> obj_ptr
224  = ROL::makePtr<NonlinearLeastSquaresObjective_Dynamic<Real>>(con_ptr,c,uo_ptr,z_ptr,ts_ptr,true);
225  ROL::ParameterList parlist;
226  parlist.sublist("Status Test").set("Gradient Tolerance",ctol);
227  parlist.sublist("Status Test").set("Step Tolerance",stol_);
228  parlist.sublist("Status Test").set("Iteration Limit",maxit_);
229  parlist.sublist("Step").sublist("Trust Region").set("Subproblem Solver","Truncated CG");
230  parlist.sublist("General").sublist("Krylov").set("Iteration Limit",100);
231  ROL::Ptr<Algorithm<Real>> algo = ROL::makePtr<Algorithm<Real>>("Trust Region",parlist,false);
232  algo->run(un,*obj_ptr,print_);
233  value(c,uo,un,z,ts);
234  }
235  if (solverType_==2 || (solverType_==4 && cnorm > ctol)) {
236  ROL::Ptr<DynamicConstraint<Real>> con_ptr = ROL::makePtrFromRef(*this);
237  ROL::Ptr<const Vector<Real>> uo_ptr = ROL::makePtrFromRef(uo);
238  ROL::Ptr<const Vector<Real>> z_ptr = ROL::makePtrFromRef(z);
239  ROL::Ptr<const TimeStamp<Real>> ts_ptr = ROL::makePtrFromRef(ts);
240  ROL::Ptr<Constraint<Real>> cn_ptr
241  = ROL::makePtr<Constraint_DynamicState<Real>>(con_ptr,uo_ptr,z_ptr,ts_ptr);
242  ROL::Ptr<Objective<Real>> obj_ptr
243  = ROL::makePtr<Objective_FSsolver<Real>>();
244  ROL::ParameterList parlist;
245  parlist.sublist("Status Test").set("Constraint Tolerance",ctol);
246  parlist.sublist("Status Test").set("Step Tolerance",stol_);
247  parlist.sublist("Status Test").set("Iteration Limit",maxit_);
248  ROL::Ptr<Algorithm<Real>> algo = ROL::makePtr<Algorithm<Real>>("Composite Step",parlist,false);
249  ROL::Ptr<Vector<Real>> l = c.dual().clone();
250  algo->run(un,*l,*obj_ptr,*cn_ptr,print_);
251  value(c,uo,un,z,ts);
252  }
253  if (solverType_ > 4 || solverType_ < 0) {
254  ROL_TEST_FOR_EXCEPTION(true, std::invalid_argument,
255  ">>> ERROR (ROL:DynamicConstraint:solve): Invalid solver type!");
256  }
257  }
258 
279  virtual void setSolveParameters(ROL::ParameterList &parlist) {
280  ROL::ParameterList & list = parlist.sublist("Dynamic Constraint").sublist("Solve");
281  atol_ = list.get("Absolute Residual Tolerance", DEFAULT_atol_);
282  rtol_ = list.get("Relative Residual Tolerance", DEFAULT_rtol_);
283  maxit_ = list.get("Iteration Limit", DEFAULT_maxit_);
284  decr_ = list.get("Sufficient Decrease Tolerance", DEFAULT_decr_);
285  stol_ = list.get("Step Tolerance", DEFAULT_stol_);
286  factor_ = list.get("Backtracking Factor", DEFAULT_factor_);
287  print_ = list.get("Output Iteration History", DEFAULT_print_);
288  zero_ = list.get("Zero Initial Guess", DEFAULT_zero_);
289  solverType_ = list.get("Solver Type", DEFAULT_solverType_);
290  }
291 
292  //----------------------------------------------------------------------------
293  // Partial Jacobians
294  virtual void applyJacobian_uo( V& jv, const V& vo, const V& uo,
295  const V& un, const V& z,
296  const TS& ts ) const {}
297 
298  virtual void applyJacobian_un( V& jv, const V& vn, const V& uo,
299  const V& un, const V& z,
300  const TS& ts ) const {}
301 
302  virtual void applyJacobian_z( V& jv, const V& vz, const V& uo,
303  const V& un, const V& z,
304  const TS& ts ) const {}
305 
306  //----------------------------------------------------------------------------
307  // Adjoint partial Jacobians
308  virtual void applyAdjointJacobian_uo( V& ajv, const V& dualv, const V& uo,
309  const V& un, const V& z,
310  const TS& ts ) const {}
311 
312  virtual void applyAdjointJacobian_un( V& ajv, const V& dualv, const V& uo,
313  const V& un, const V& z,
314  const TS& ts ) const {}
315 
316  virtual void applyAdjointJacobian_z( V& ajv, const V& dualv, const V& uo,
317  const V& un, const V& z,
318  const TS& ts ) const {}
319 
320  //----------------------------------------------------------------------------
321  // Inverses
322  virtual void applyInverseJacobian_un( V& ijv, const V& vn, const V& uo,
323  const V& un, const V& z,
324  const TS& ts ) const {}
325 
326  virtual void applyInverseAdjointJacobian_un( V& iajv, const V& vn, const V& uo,
327  const V& un, const V& z,
328  const TS& ts ) const {}
329 
330  //----------------------------------------------------------------------------
331  // Adjoint Hessian components
332  virtual void applyAdjointHessian_un_un( V& ahwv, const V& wn, const V& vn,
333  const V& uo, const V& un,
334  const V& z, const TS& ts ) const {
335  ahwv.zero();
336  }
337 
338  virtual void applyAdjointHessian_un_uo( V& ahwv, const V& w, const V& vn,
339  const V& uo, const V& un,
340  const V& z, const TS& ts ) const {
341  ahwv.zero();
342  }
343 
344  virtual void applyAdjointHessian_un_z( V& ahwv, const V& w, const V& vn,
345  const V& uo, const V& un,
346  const V& z, const TS& ts ) const {
347  ahwv.zero();
348  }
349 
350  virtual void applyAdjointHessian_uo_un( V& ahwv, const V& w, const V& vo,
351  const V& uo, const V& un,
352  const V& z, const TS& ts ) const {
353  ahwv.zero();
354  }
355 
356  virtual void applyAdjointHessian_uo_uo( V& ahwv, const V& w, const V& v,
357  const V& uo, const V& un,
358  const V& z, const TS& ts ) const {
359  ahwv.zero();
360  }
361 
362  virtual void applyAdjointHessian_uo_z( V& ahwv, const V& w, const V& vo,
363  const V& uo, const V& un,
364  const V& z, const TS& ts ) const {
365  ahwv.zero();
366  }
367 
368  virtual void applyAdjointHessian_z_un( V& ahwv, const V& w, const V& vz,
369  const V& uo, const V& un,
370  const V& z, const TS& ts ) const {
371  ahwv.zero();
372  }
373 
374  virtual void applyAdjointHessian_z_uo( V& ahwv, const V& w, const V& vz,
375  const V& uo, const V& un,
376  const V& z, const TS& ts ) const {
377  ahwv.zero();
378  }
379 
380  virtual void applyAdjointHessian_z_z( V& ahwv, const V& w, const V& vz,
381  const V& uo, const V& un,
382  const V& z, const TS& ts ) const {
383  ahwv.zero();
384  }
385 
386 }; // DynamicConstraint
387 
388 
389 } // namespace ROL
390 
391 
392 #endif // ROL_DYNAMICCONSTRAINT_HPP
393 
Ptr< Vector< Real > > jv_
virtual void update_un(const V &un, const TS &ts)
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:226
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
Contains local time step information.
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
virtual void setSolveParameters(ROL::ParameterList &parlist)
Set solve parameters.
virtual void applyAdjointJacobian_z(V &ajv, const V &dualv, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void applyInverseJacobian_un(V &ijv, const V &vn, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void value(V &c, const V &uo, const V &un, const V &z, const TS &ts) const =0
virtual void applyAdjointJacobian_uo(V &ajv, const V &dualv, const V &uo, const V &un, const V &z, const TS &ts) const
virtual void update_z(const V &z, const TS &ts)
virtual void update_uo(const V &uo, const TS &ts)
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:209
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
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:91