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
50
54#include "ROL_Algorithm.hpp"
55#include "ROL_CompositeStep.hpp"
58#include "ROL_Types.hpp"
59
78namespace ROL {
79
80template<typename Real>
81class DynamicConstraint : public DynamicFunction<Real> {
82private:
83 // Additional vector storage for solve
84 Ptr<Vector<Real>> unew_;
85 Ptr<Vector<Real>> jv_;
86
87 // Default parameters for solve (backtracking Newton)
88 const Real DEFAULT_atol_;
89 const Real DEFAULT_rtol_;
90 const Real DEFAULT_stol_;
91 const Real DEFAULT_factor_;
92 const Real DEFAULT_decr_;
93 const int DEFAULT_maxit_;
94 const bool DEFAULT_print_;
95 const bool DEFAULT_zero_;
97
98 // User-set parameters for solve (backtracking Newton)
99 Real atol_;
100 Real rtol_;
101 Real stol_;
103 Real decr_;
105 bool print_;
106 bool zero_;
108
109 // Flag to initialize vector storage in solve
111
112public:
113
117
119
120 DynamicConstraint( std::initializer_list<std::string> zero_deriv_terms={} ):
121 DynamicFunction<Real>(zero_deriv_terms),
122 unew_ ( nullPtr ),
123 jv_ ( nullPtr ),
124 DEFAULT_atol_ ( 1e-4*std::sqrt(ROL_EPSILON<Real>()) ),
125 DEFAULT_rtol_ ( 1e0 ),
126 DEFAULT_stol_ ( std::sqrt(ROL_EPSILON<Real>()) ),
127 DEFAULT_factor_ ( 0.5 ),
128 DEFAULT_decr_ ( 1e-4 ),
129 DEFAULT_maxit_ ( 500 ),
130 DEFAULT_print_ ( false ),
131 DEFAULT_zero_ ( false ),
142 firstSolve_ ( true ) {}
143
144
145
146 virtual void update( const V& uo, const V& un, const V& z, const TS& ts ) {
147 update_uo( uo, ts );
148 update_un( un, ts );
149 update_z( z, ts );
150 }
151
152 using DynamicFunction<Real>::update_uo;
153 using DynamicFunction<Real>::update_un;
154 using DynamicFunction<Real>::update_z;
155
156 virtual void value( V& c, const V& uo, const V& un,
157 const V& z, const TS& ts ) const = 0;
158
159 virtual void solve( V& c, const V& uo, V& un,
160 const V& z, const TS& ts ) {
161 if ( zero_ ) {
162 un.zero();
163 }
164 update(uo,un,z,ts);
165 value(c,uo,un,z,ts);
166 Real cnorm = c.norm();
167 const Real ctol = std::min(atol_, rtol_*cnorm);
168 if (solverType_==0 || solverType_==3 || solverType_==4) {
169 if ( firstSolve_ ) {
170 unew_ = un.clone();
171 jv_ = un.clone();
172 firstSolve_ = false;
173 }
174 Real alpha(1), tmp(0);
175 int cnt = 0;
176 if ( print_ ) {
177 std::cout << "\n Default DynamicConstraint::solve\n";
178 std::cout << " ";
179 std::cout << std::setw(6) << std::left << "iter";
180 std::cout << std::setw(15) << std::left << "rnorm";
181 std::cout << std::setw(15) << std::left << "alpha";
182 std::cout << "\n";
183 }
184 for (cnt = 0; cnt < maxit_; ++cnt) {
185 // Compute Newton step
186 applyInverseJacobian_un(*jv_,c,uo,un,z,ts);
187 unew_->set(un);
188 unew_->axpy(-alpha, *jv_);
189 //update_un(*unew_,ts);
190 update(uo,*unew_,z,ts);
191 value(c,uo,*unew_,z,ts);
192 tmp = c.norm();
193 // Perform backtracking line search
194 while ( tmp > (1.0-decr_*alpha)*cnorm &&
195 alpha > stol_ ) {
196 alpha *= factor_;
197 unew_->set(un);
198 unew_->axpy(-alpha,*jv_);
199 //update_un(*unew_,ts);
200 update(uo,*unew_,z,ts);
201 value(c,uo,*unew_,z,ts);
202 tmp = c.norm();
203 }
204 if ( print_ ) {
205 std::cout << " ";
206 std::cout << std::setw(6) << std::left << cnt;
207 std::cout << std::scientific << std::setprecision(6);
208 std::cout << std::setw(15) << std::left << tmp;
209 std::cout << std::scientific << std::setprecision(6);
210 std::cout << std::setw(15) << std::left << alpha;
211 std::cout << "\n";
212 }
213 // Update iterate
214 cnorm = tmp;
215 un.set(*unew_);
216 if (cnorm <= ctol) { // = covers the case of identically zero residual
217 break;
218 }
219 update(uo,un,z,ts);
220 alpha = 1.0;
221 }
222 }
223 if (solverType_==1 || (solverType_==3 && cnorm > ctol)) {
224 Ptr<DynamicConstraint<Real>> con_ptr = makePtrFromRef(*this);
225 Ptr<const Vector<Real>> uo_ptr = makePtrFromRef(uo);
226 Ptr<const Vector<Real>> z_ptr = makePtrFromRef(z);
227 Ptr<const TimeStamp<Real>> ts_ptr = makePtrFromRef(ts);
228 Ptr<Objective<Real>> obj_ptr
229 = makePtr<NonlinearLeastSquaresObjective_Dynamic<Real>>(con_ptr,c,uo_ptr,z_ptr,ts_ptr,true);
230 ParameterList parlist;
231 parlist.sublist("Status Test").set("Gradient Tolerance",ctol);
232 parlist.sublist("Status Test").set("Step Tolerance",stol_);
233 parlist.sublist("Status Test").set("Iteration Limit",maxit_);
234 parlist.sublist("Step").sublist("Trust Region").set("Subproblem Solver","Truncated CG");
235 parlist.sublist("General").sublist("Krylov").set("Iteration Limit",100);
236 Ptr<Step<Real>> step = makePtr<TrustRegionStep<Real>>(parlist);
237 Ptr<StatusTest<Real>> status = makePtr<StatusTest<Real>>(parlist);
238 Ptr<Algorithm<Real>> algo = makePtr<Algorithm<Real>>(step,status,false);
239 algo->run(un,*obj_ptr,print_);
240 value(c,uo,un,z,ts);
241 }
242 if (solverType_==2 || (solverType_==4 && cnorm > ctol)) {
243 Ptr<DynamicConstraint<Real>> con_ptr = makePtrFromRef(*this);
244 Ptr<const Vector<Real>> uo_ptr = makePtrFromRef(uo);
245 Ptr<const Vector<Real>> z_ptr = makePtrFromRef(z);
246 Ptr<const TimeStamp<Real>> ts_ptr = makePtrFromRef(ts);
247 Ptr<Constraint<Real>> cn_ptr
248 = makePtr<Constraint_DynamicState<Real>>(con_ptr,uo_ptr,z_ptr,ts_ptr);
249 Ptr<Objective<Real>> obj_ptr
250 = makePtr<Objective_FSsolver<Real>>();
251 ParameterList parlist;
252 parlist.sublist("Status Test").set("Constraint Tolerance",ctol);
253 parlist.sublist("Status Test").set("Step Tolerance",stol_);
254 parlist.sublist("Status Test").set("Iteration Limit",maxit_);
255 Ptr<Step<Real>> step = makePtr<CompositeStep<Real>>(parlist);
256 Ptr<StatusTest<Real>> status = makePtr<ConstraintStatusTest<Real>>(parlist);
257 Ptr<Algorithm<Real>> algo = makePtr<Algorithm<Real>>(step,status,false);
258 Ptr<Vector<Real>> l = c.dual().clone();
259 algo->run(un,*l,*obj_ptr,*cn_ptr,print_);
260 value(c,uo,un,z,ts);
261 }
262 if (solverType_ > 4 || solverType_ < 0) {
263 ROL_TEST_FOR_EXCEPTION(true, std::invalid_argument,
264 ">>> ERROR (ROL:DynamicConstraint:solve): Invalid solver type!");
265 }
266 }
267
288 virtual void setSolveParameters(ParameterList &parlist) {
289 ParameterList & list = parlist.sublist("Dynamic Constraint").sublist("Solve");
290 atol_ = list.get("Absolute Residual Tolerance", DEFAULT_atol_);
291 rtol_ = list.get("Relative Residual Tolerance", DEFAULT_rtol_);
292 maxit_ = list.get("Iteration Limit", DEFAULT_maxit_);
293 decr_ = list.get("Sufficient Decrease Tolerance", DEFAULT_decr_);
294 stol_ = list.get("Step Tolerance", DEFAULT_stol_);
295 factor_ = list.get("Backtracking Factor", DEFAULT_factor_);
296 print_ = list.get("Output Iteration History", DEFAULT_print_);
297 zero_ = list.get("Zero Initial Guess", DEFAULT_zero_);
298 solverType_ = list.get("Solver Type", DEFAULT_solverType_);
299 }
300
301 //----------------------------------------------------------------------------
302 // Partial Jacobians
303 virtual void applyJacobian_uo( V& jv, const V& vo, const V& uo,
304 const V& un, const V& z,
305 const TS& ts ) const {}
306
307 virtual void applyJacobian_un( V& jv, const V& vn, const V& uo,
308 const V& un, const V& z,
309 const TS& ts ) const {}
310
311 virtual void applyJacobian_z( V& jv, const V& vz, const V& uo,
312 const V& un, const V& z,
313 const TS& ts ) const {}
314
315 //----------------------------------------------------------------------------
316 // Adjoint partial Jacobians
317 virtual void applyAdjointJacobian_uo( V& ajv, const V& dualv, const V& uo,
318 const V& un, const V& z,
319 const TS& ts ) const {}
320
321 virtual void applyAdjointJacobian_un( V& ajv, const V& dualv, const V& uo,
322 const V& un, const V& z,
323 const TS& ts ) const {}
324
325 virtual void applyAdjointJacobian_z( V& ajv, const V& dualv, const V& uo,
326 const V& un, const V& z,
327 const TS& ts ) const {}
328
329 //----------------------------------------------------------------------------
330 // Inverses
331 virtual void applyInverseJacobian_un( V& ijv, const V& vn, const V& uo,
332 const V& un, const V& z,
333 const TS& ts ) const {}
334
335 virtual void applyInverseAdjointJacobian_un( V& iajv, const V& vn, const V& uo,
336 const V& un, const V& z,
337 const TS& ts ) const {}
338
339 //----------------------------------------------------------------------------
340 // Adjoint Hessian components
341 virtual void applyAdjointHessian_un_un( V& ahwv, const V& wn, const V& vn,
342 const V& uo, const V& un,
343 const V& z, const TS& ts ) const {
344 ahwv.zero();
345 }
346
347 virtual void applyAdjointHessian_un_uo( V& ahwv, const V& w, const V& vn,
348 const V& uo, const V& un,
349 const V& z, const TS& ts ) const {
350 ahwv.zero();
351 }
352
353 virtual void applyAdjointHessian_un_z( V& ahwv, const V& w, const V& vn,
354 const V& uo, const V& un,
355 const V& z, const TS& ts ) const {
356 ahwv.zero();
357 }
358
359 virtual void applyAdjointHessian_uo_un( V& ahwv, const V& w, const V& vo,
360 const V& uo, const V& un,
361 const V& z, const TS& ts ) const {
362 ahwv.zero();
363 }
364
365 virtual void applyAdjointHessian_uo_uo( V& ahwv, const V& w, const V& v,
366 const V& uo, const V& un,
367 const V& z, const TS& ts ) const {
368 ahwv.zero();
369 }
370
371 virtual void applyAdjointHessian_uo_z( V& ahwv, const V& w, const V& vo,
372 const V& uo, const V& un,
373 const V& z, const TS& ts ) const {
374 ahwv.zero();
375 }
376
377 virtual void applyAdjointHessian_z_un( V& ahwv, const V& w, const V& vz,
378 const V& uo, const V& un,
379 const V& z, const TS& ts ) const {
380 ahwv.zero();
381 }
382
383 virtual void applyAdjointHessian_z_uo( V& ahwv, const V& w, const V& vz,
384 const V& uo, const V& un,
385 const V& z, const TS& ts ) const {
386 ahwv.zero();
387 }
388
389 virtual void applyAdjointHessian_z_z( V& ahwv, const V& w, const V& vz,
390 const V& uo, const V& un,
391 const V& z, const TS& ts ) const {
392 ahwv.zero();
393 }
394
395}; // DynamicConstraint
396
397
398} // namespace ROL
399
400
401#endif // ROL_DYNAMICCONSTRAINT_HPP
402
Contains definitions of custom data types in ROL.
Defines the time-dependent constraint operator interface for simulation-based optimization.
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 setSolveParameters(ParameterList &parlist)
Set solve parameters.
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_uo_un(V &ahwv, const V &w, const V &vo, const V &uo, const V &un, const V &z, const TS &ts) const
Ptr< Vector< Real > > unew_
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 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 applyInverseAdjointJacobian_un(V &iajv, 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 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 void solve(V &c, const V &uo, V &un, const V &z, const TS &ts)
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 value(V &c, const V &uo, const V &un, const V &z, const TS &ts) const =0
virtual void applyJacobian_uo(V &jv, const V &vo, const V &uo, const V &un, const V &z, const TS &ts) const
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(const V &uo, const V &un, const V &z, const TS &ts)
virtual void applyInverseJacobian_un(V &ijv, const V &vn, const V &uo, const V &un, const V &z, const TS &ts) const
Ptr< Vector< Real > > jv_
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
virtual void applyJacobian_un(V &jv, const V &vn, 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 applyAdjointJacobian_z(V &ajv, const V &dualv, const V &uo, const V &un, const V &z, const TS &ts) const
DynamicConstraint(std::initializer_list< std::string > zero_deriv_terms={})
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
Provides update interface, casting and vector management to DynamicConstraint and DynamicObjective.
virtual void update_uo(const V &x, const TS &ts)
virtual void update_z(const V &x, const TS &ts)
virtual void update_un(const V &x, const TS &ts)
Defines the linear algebra of vector space on a generic partitioned vector.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual Real norm() const =0
Returns where .
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
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 zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Contains local time step information.