ROL
ROL_ReducedDynamicObjective.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#ifndef ROL_REDUCEDDYNAMICOBJECTIVE_HPP
45#define ROL_REDUCEDDYNAMICOBJECTIVE_HPP
46
47#include "ROL_Ptr.hpp"
48#include "ROL_Sketch.hpp"
49#include "ROL_Objective.hpp"
52
53#include <fstream>
54
77namespace ROL {
78
79template<typename Real>
80class ReducedDynamicObjective : public Objective<Real> {
81 using size_type = typename std::vector<Real>::size_type;
82private:
83 // Problem data.
84 const Ptr<DynamicObjective<Real>> obj_;
85 const Ptr<DynamicConstraint<Real>> con_;
86 const Ptr<Vector<Real>> u0_;
87 const std::vector<TimeStamp<Real>> timeStamp_;
89 // General sketch information.
90 const bool useSketch_;
91 // State sketch information.
93 Ptr<Sketch<Real>> stateSketch_;
94 // Adjoint sketch information.
96 Ptr<Sketch<Real>> adjointSketch_;
97 // State sensitivity sketch information.
99 Ptr<Sketch<Real>> stateSensSketch_;
100 // Inexactness information.
101 const bool useInexact_;
104 const bool syncHessRank_;
105 // Vector storage for intermediate computations.
106 std::vector<Ptr<Vector<Real>>> uhist_; // State history.
107 std::vector<Ptr<Vector<Real>>> lhist_; // Adjoint history.
108 std::vector<Ptr<Vector<Real>>> whist_; // State sensitivity history.
109 std::vector<Ptr<Vector<Real>>> phist_; // Adjoint sensitivity history.
110 Ptr<Vector<Real>> cprimal_; // Primal constraint vector.
111 Ptr<Vector<Real>> crhs_; // Primal constraint vector.
112 Ptr<Vector<Real>> udual_; // Dual state vector.
113 Ptr<Vector<Real>> rhs_; // Dual state vector.
114 Ptr<Vector<Real>> zdual_; // Dual control vector.
115 Real val_; // Objective function value.
116 // Flags.
117 bool isValueComputed_; // Whether objective has been evaluated.
118 bool isStateComputed_; // Whether state has been solved.
119 bool isAdjointComputed_; // Whether adjoint has been solved.
120 bool useHessian_; // Whether to use Hessian or FD.
121 bool useSymHess_; // Whether to use symmetric Hessian approximation.
122 // Output information
123 Ptr<std::ostream> stream_;
124 const bool print_;
125 const int freq_;
126
128 return static_cast<PartitionedVector<Real>&>(x);
129 }
130
132 return static_cast<const PartitionedVector<Real>&>(x);
133 }
134
135 void throwError(const int err, const std::string &sfunc,
136 const std::string &func, const int line) const {
137 if (err != 0) {
138 std::stringstream out;
139 out << ">>> ROL::ReducedDynamicObjective::" << func << " (Line " << line << "): ";
140 if (err == 1) {
141 out << "Reconstruction has already been called!";
142 }
143 else if (err == 2) {
144 out << "Input column index exceeds total number of columns!";
145 }
146 else if (err == 3) {
147 out << "Reconstruction failed to compute domain-space basis!";
148 }
149 else if (err == 4) {
150 out << "Reconstruction failed to compute range-space basis!";
151 }
152 else if (err == 5) {
153 out << "Reconstruction failed to generate core matrix!";
154 }
155 throw Exception::NotImplemented(out.str());
156 }
157 }
158
159public:
161 const Ptr<DynamicConstraint<Real>> &con,
162 const Ptr<Vector<Real>> &u0,
163 const Ptr<Vector<Real>> &zvec,
164 const Ptr<Vector<Real>> &cvec,
165 const std::vector<TimeStamp<Real>> &timeStamp,
166 ROL::ParameterList &pl,
167 const Ptr<std::ostream> &stream = nullPtr)
168 : obj_ ( obj ), // Dynamic objective function.
169 con_ ( con ), // Dynamic constraint function.
170 u0_ ( u0 ), // Initial condition.
171 timeStamp_ ( timeStamp ), // Vector of time stamps.
172 Nt_ ( timeStamp.size() ), // Number of time intervals.
173 useSketch_ ( pl.get("Use Sketching", false) ), // Use state sketch if true.
174 rankState_ ( pl.get("State Rank", 10) ), // Rank of state sketch.
175 stateSketch_ ( nullPtr ), // State sketch object.
176 rankAdjoint_ ( pl.get("Adjoint Rank", 10) ), // Rank of adjoint sketch.
177 adjointSketch_ ( nullPtr ), // Adjoint sketch object.
178 rankStateSens_ ( pl.get("State Sensitivity Rank", 10) ), // Rank of state sensitivity sketch.
179 stateSensSketch_ ( nullPtr ), // State sensitivity sketch object.
180 useInexact_ ( pl.get("Adaptive Rank", false) ), // Update rank adaptively.
181 updateFactor_ ( pl.get("Rank Update Factor", 2) ), // Rank update factor.
182 maxRank_ ( pl.get("Maximum Rank", 100) ), // Maximum rank.
183 syncHessRank_ ( pl.get("Sync Hessian Rank", useInexact_) ), // Sync rank for Hessian storage.
184 isValueComputed_ ( false ), // Flag indicating whether value has been computed.
185 isStateComputed_ ( false ), // Flag indicating whether state has been computed.
186 isAdjointComputed_ ( false ), // Flag indicating whether adjoint has been computed.
187 useHessian_ ( pl.get("Use Hessian", true) ), // Flag indicating whether to use the Hessian.
188 useSymHess_ ( pl.get("Use Only Sketched Sensitivity", true) ), // Flag indicating whether to use symmetric sketched Hessian.
189 stream_ ( stream ), // Output stream to print sketch information.
190 print_ ( pl.get("Print Optimization Vector", false) ), // Print control vector to file
191 freq_ ( pl.get("Output Frequency", 5) ) { // Print frequency for control vector
192 uhist_.clear(); lhist_.clear(); whist_.clear(); phist_.clear();
193 if (useSketch_) { // Only maintain a sketch of the state time history
194 Real orthTol = pl.get("Orthogonality Tolerance", 1e2*ROL_EPSILON<Real>());
195 int orthIt = pl.get("Reorthogonalization Iterations", 5);
196 bool trunc = pl.get("Truncate Approximation", false);
197 if (syncHessRank_) {
200 }
201 unsigned dseed = pl.get("State Domain Seed",0);
202 unsigned rseed = pl.get("State Range Seed",0);
203 stateSketch_ = makePtr<Sketch<Real>>(*u0_,static_cast<int>(Nt_)-1,rankState_,
204 orthTol,orthIt,trunc,dseed,rseed);
205 uhist_.push_back(u0_->clone());
206 uhist_.push_back(u0_->clone());
207 lhist_.push_back(cvec->dual().clone());
208 dseed = pl.get("Adjoint Domain Seed",0);
209 rseed = pl.get("Adjoint Range Seed",0);
210 adjointSketch_ = makePtr<Sketch<Real>>(*u0_,static_cast<int>(Nt_)-1,rankAdjoint_,
211 orthTol,orthIt,trunc,dseed,rseed);
212 if (useHessian_) {
213 dseed = pl.get("State Sensitivity Domain Seed",0);
214 rseed = pl.get("State Sensitivity Range Seed",0);
215 stateSensSketch_ = makePtr<Sketch<Real>>(*u0_,static_cast<int>(Nt_)-1,
216 rankStateSens_,orthTol,orthIt,trunc,dseed,rseed);
217 whist_.push_back(u0_->clone());
218 whist_.push_back(u0_->clone());
219 phist_.push_back(cvec->dual().clone());
220 }
221 }
222 else { // Store entire state time history
223 for (size_type i = 0; i < Nt_; ++i) {
224 uhist_.push_back(u0_->clone());
225 lhist_.push_back(cvec->dual().clone());
226 if (useHessian_) {
227 whist_.push_back(u0_->clone());
228 phist_.push_back(cvec->dual().clone());
229 }
230 }
231 }
232 cprimal_ = cvec->clone();
233 crhs_ = cvec->clone();
234 udual_ = u0_->dual().clone();
235 rhs_ = u0_->dual().clone();
236 zdual_ = zvec->dual().clone();
237 }
238
239 Ptr<Vector<Real>> makeDynamicVector(const Vector<Real> &x) const {
241 }
242
243 void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
244 if (useSketch_) {
245 stateSketch_->update();
246 if (flag == true) {
247 adjointSketch_->update();
248 if (useHessian_) {
249 stateSensSketch_->update();
250 }
251 }
252 }
253 for (size_type i = 0; i < uhist_.size(); ++i) {
254 uhist_[i]->zero();
255 }
256 val_ = static_cast<Real>(0);
257 isValueComputed_ = false;
258 isStateComputed_ = false;
259 if (flag == true) {
260 isAdjointComputed_ = false;
261 }
262 if (iter >= 0 && iter%freq_==0 && print_) {
263 std::stringstream name;
264 name << "optvector." << iter << ".txt";
265 std::ofstream file;
266 file.open(name.str());
267 x.print(file);
268 file.close();
269 }
270 }
271
272 Real value( const Vector<Real> &x, Real &tol ) {
273 if (!isValueComputed_) {
274 int eflag(0);
275 const Real one(1);
276 const PartitionedVector<Real> &xp = partition(x);
277 // Set initial condition
278 uhist_[0]->set(*u0_);
279 if (useSketch_) {
280 stateSketch_->update();
281 uhist_[1]->set(*u0_);
282 }
283 // Run time stepper
284 Real valk(0);
285 size_type index;
286 for (size_type k = 1; k < Nt_; ++k) {
287 index = (useSketch_ ? 1 : k);
288 if (!useSketch_) {
289 uhist_[index]->set(*uhist_[index-1]);
290 }
291 // Update dynamic constraint
292 con_->update_uo(*uhist_[index-1], timeStamp_[k]);
293 con_->update_z(*xp.get(k), timeStamp_[k]);
294 // Solve state on current time interval
295 con_->solve(*cprimal_, *uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
296 // Update dynamic objective
297 obj_->update(*uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
298 // Compute objective function value on current time interval
299 valk = obj_->value(*uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
300 // Update total objective function value
301 val_ += valk;
302 // Sketch state
303 if (useSketch_) {
304 eflag = stateSketch_->advance(one,*uhist_[1],static_cast<int>(k)-1,one);
305 throwError(eflag,"advance","value",273);
306 uhist_[0]->set(*uhist_[1]);
307 }
308 }
309 isValueComputed_ = true;
310 isStateComputed_ = true;
311 }
312 return val_;
313 }
314
315 void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
316 int eflag(0);
318 gp.get(0)->zero(); // zero for the nonexistant zeroth control interval
319 const PartitionedVector<Real> &xp = partition(x);
320 const Real one(1);
321 size_type uindex = (useSketch_ ? 1 : Nt_-1);
322 size_type lindex = (useSketch_ ? 0 : Nt_-1);
323 // Must first compute the value
324 solveState(x);
325 if (useSketch_ && useInexact_) {
326 if (stream_ != nullPtr) {
327 *stream_ << std::string(80,'=') << std::endl;
328 *stream_ << " ROL::ReducedDynamicObjective::gradient" << std::endl;
329 }
330 tol = updateSketch(x,tol);
331 if (stream_ != nullPtr) {
332 *stream_ << " State Rank for Gradient Computation: " << rankState_ << std::endl;
333 *stream_ << " Residual Norm: " << tol << std::endl;
334 *stream_ << std::string(80,'=') << std::endl;
335 }
336 }
337 // Recover state from sketch
338 if (useSketch_) {
339 uhist_[1]->set(*uhist_[0]);
340 eflag = stateSketch_->reconstruct(*uhist_[0],static_cast<int>(Nt_)-3);
341 throwError(eflag,"reconstruct","gradient",309);
342 if (isAdjointComputed_) {
343 eflag = adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(Nt_)-2);
344 throwError(eflag,"reconstruct","gradient",312);
345 }
346 }
347 // Update dynamic constraint and objective
348 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
349 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
350 // Solve for terminal condition
351 if (!isAdjointComputed_) {
353 *uhist_[uindex-1], *uhist_[uindex],
354 *xp.get(Nt_-1), timeStamp_[Nt_-1]);
355 if (useSketch_) {
356 eflag = adjointSketch_->advance(one,*lhist_[0],static_cast<int>(Nt_)-2,one);
357 throwError(eflag,"advance","gradient",325);
358 }
359 }
360 // Update gradient on terminal interval
361 updateGradient(*gp.get(Nt_-1), *lhist_[lindex],
362 *uhist_[uindex-1], *uhist_[uindex],
363 *xp.get(Nt_-1), timeStamp_[Nt_-1]);
364 // Run reverse time stepper
365 for (size_type k = Nt_-2; k > 0; --k) {
366 if (!isAdjointComputed_) {
367 // Compute k+1 component of rhs
368 computeAdjointRHS(*rhs_, *lhist_[lindex],
369 *uhist_[uindex-1], *uhist_[uindex],
370 *xp.get(k+1), timeStamp_[k+1]);
371 }
372 uindex = (useSketch_ ? 1 : k);
373 lindex = (useSketch_ ? 0 : k);
374 // Recover state from sketch
375 if (useSketch_) {
376 uhist_[1]->set(*uhist_[0]);
377 if (k==1) {
378 uhist_[0]->set(*u0_);
379 }
380 else {
381 eflag = stateSketch_->reconstruct(*uhist_[0],static_cast<int>(k)-2);
382 throwError(eflag,"reconstruct","gradient",350);
383 }
384 if (isAdjointComputed_) {
385 eflag = adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(k)-1);
386 throwError(eflag,"reconstruct","gradient",354);
387 }
388 }
389 // Update dynamic constraint and objective
390 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
391 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
392 // Solve for adjoint on interval k
393 if (!isAdjointComputed_) {
394 advanceAdjoint(*lhist_[lindex], *rhs_,
395 *uhist_[uindex-1], *uhist_[uindex],
396 *xp.get(k), timeStamp_[k]);
397 if (useSketch_) {
398 eflag = adjointSketch_->advance(one,*lhist_[0],static_cast<int>(k)-1,one);
399 throwError(eflag,"advance","gradient",367);
400 }
401 }
402 // Update gradient on interval k
403 updateGradient(*gp.get(k), *lhist_[lindex],
404 *uhist_[uindex-1], *uhist_[uindex],
405 *xp.get(k), timeStamp_[k]);
406 }
407 isAdjointComputed_ = true;
408 }
409
410 void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
411 int eflag(0);
412 if (useHessian_) {
413 // Must first solve the state and adjoint equations
414 solveState(x);
415 solveAdjoint(x);
416 // Now compute Hessian
417 const Real one(1);
418 const PartitionedVector<Real> &xp = partition(x);
419 const PartitionedVector<Real> &vp = partition(v);
421 hvp.get(0)->zero(); // zero for the nonexistant zeroth control interval
422 // Compute state sensitivity
423 whist_[0]->zero();
424 if (useSketch_) {
425 stateSensSketch_->update();
426 //stateSensSketch_->advance(one,*whist_[0],0,one);
427 uhist_[0]->set(*u0_);
428 }
429 size_type uindex, lindex;
430 for (size_type k = 1; k < Nt_; ++k) {
431 uindex = (useSketch_ ? 1 : k);
432 lindex = (useSketch_ ? 0 : k);
433 // Reconstruct sketched state
434 if (useSketch_) {
435 eflag = stateSketch_->reconstruct(*uhist_[1],static_cast<int>(k)-1);
436 throwError(eflag,"reconstruct","hessVec",404);
437 eflag = adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(k)-1);
438 throwError(eflag,"reconstruct","hessVec",406);
439 }
440 // Update dynamic constraint and objective
441 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
442 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
443 // Advance state sensitivity on current time interval
444 advanceStateSens(*whist_[uindex],
445 *vp.get(k), *whist_[uindex-1],
446 *uhist_[uindex-1], *uhist_[uindex],
447 *xp.get(k), timeStamp_[k]);
448 // Set control only Hessian
450 *vp.get(k), *lhist_[lindex],
451 *uhist_[uindex-1], *uhist_[uindex],
452 *xp.get(k), timeStamp_[k]);
453 if (!useSymHess_) {
454 // Add mixed derivative Hessian
455 addMixedHessLag(*hvp.get(k), *lhist_[lindex],
456 *whist_[uindex-1], *whist_[uindex],
457 *uhist_[uindex-1], *uhist_[uindex],
458 *xp.get(k), timeStamp_[k]);
459 }
460 // Sketch state sensitivity
461 if (useSketch_) {
462 eflag = stateSensSketch_->advance(one,*whist_[1],static_cast<int>(k)-1,one);
463 throwError(eflag,"advance","hessVec",431);
464 whist_[0]->set(*whist_[1]);
465 uhist_[0]->set(*uhist_[1]);
466 }
467 }
468
469 // Compute adjoint sensitivity
470 uindex = (useSketch_ ? 1 : Nt_-1);
471 lindex = (useSketch_ ? 0 : Nt_-1);
472 if (useSketch_) {
473 // Recover terminal state
474 uhist_[1]->set(*uhist_[0]);
475 eflag = stateSketch_->reconstruct(*uhist_[0],static_cast<int>(Nt_)-3);
476 throwError(eflag,"reconstruct","hessVec",444);
477 // Recover terminal adjoint
478 eflag = adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(Nt_)-2);
479 throwError(eflag,"reconstruct","hessVec",447);
480 // Recover terminal state sensitivity
481 whist_[1]->set(*whist_[0]);
482 eflag = stateSensSketch_->reconstruct(*whist_[0],static_cast<int>(Nt_)-3);
483 throwError(eflag,"reconstruct","hessVec",451);
484 }
485 // Update dynamic constraint and objective
486 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
487 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
488 // Solve for terminal condition
490 *vp.get(Nt_-1), *lhist_[lindex],
491 *whist_[uindex-1], *whist_[uindex],
492 *uhist_[uindex-1], *uhist_[uindex],
493 *xp.get(Nt_-1), timeStamp_[Nt_-1]);
494 if (useSymHess_) {
495 // Add mixed derivative Hessian
496 addMixedHessLag(*hvp.get(Nt_-1), *lhist_[lindex],
497 *whist_[uindex-1], *whist_[uindex],
498 *uhist_[uindex-1], *uhist_[uindex],
499 *xp.get(Nt_-1), timeStamp_[Nt_-1]);
500 }
501 // Add adjoint sensitivity to Hessian
502 addAdjointSens(*hvp.get(Nt_-1), *phist_[lindex],
503 *uhist_[uindex-1], *uhist_[uindex],
504 *xp.get(Nt_-1), timeStamp_[Nt_-1]);
505 // Run reverse time stepper
506 for (size_type k = Nt_-2; k > 0; --k) {
507 // Compute new components of rhs
509 *whist_[uindex-1], *whist_[uindex],
510 *uhist_[uindex-1], *uhist_[uindex],
511 *xp.get(k+1), timeStamp_[k+1],
512 false);
514 *vp.get(k+1), *lhist_[lindex],
515 *uhist_[uindex-1], *uhist_[uindex],
516 *xp.get(k+1), timeStamp_[k+1],
517 true);
519 *uhist_[uindex-1], *uhist_[uindex],
520 *xp.get(k+1), timeStamp_[k+1],
521 true);
522 // Recover state, adjoint and state sensitivity from sketch
523 if (useSketch_) {
524 uhist_[1]->set(*uhist_[0]);
525 whist_[1]->set(*whist_[0]);
526 if (k==1) {
527 uhist_[0]->set(*u0_);
528 whist_[0]->zero();
529 }
530 else {
531 eflag = stateSketch_->reconstruct(*uhist_[0],static_cast<int>(k)-2);
532 throwError(eflag,"reconstruct","hessVec",500);
533 eflag = stateSensSketch_->reconstruct(*whist_[0],static_cast<int>(k)-2);
534 throwError(eflag,"reconstruct","hessVec",502);
535 }
536 eflag = adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(k)-1);
537 throwError(eflag,"reconstruct","hessVec",505);
538 }
539 uindex = (useSketch_ ? 1 : k);
540 lindex = (useSketch_ ? 0 : k);
541 // Update dynamic constraint and objective
542 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
543 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
544 // Compute old components of rhs
546 *whist_[uindex-1], *whist_[uindex],
547 *uhist_[uindex-1], *uhist_[uindex],
548 *xp.get(k), timeStamp_[k],
549 true);
551 *vp.get(k), *lhist_[lindex],
552 *uhist_[uindex-1], *uhist_[uindex],
553 *xp.get(k), timeStamp_[k],
554 true);
555 // Solve for adjoint on interval k
556 advanceAdjointSens(*phist_[lindex], *rhs_,
557 *uhist_[uindex-1], *uhist_[uindex],
558 *xp.get(k), timeStamp_[k]);
559 if (useSymHess_) {
560 // Add mixed derivative Hessian
561 addMixedHessLag(*hvp.get(k), *lhist_[lindex],
562 *whist_[uindex-1], *whist_[uindex],
563 *uhist_[uindex-1], *uhist_[uindex],
564 *xp.get(k), timeStamp_[k]);
565 }
566 // Add adjoint sensitivity to Hessian
567 addAdjointSens(*hvp.get(k), *phist_[lindex],
568 *uhist_[uindex-1], *uhist_[uindex],
569 *xp.get(k), timeStamp_[k]);
570 }
571 }
572 else {
573 Objective<Real>::hessVec(hv,v,x,tol);
574 }
575 }
576
577private:
578 /***************************************************************************/
579 /************ Method to solve state equation *******************************/
580 /***************************************************************************/
581 Real solveState(const Vector<Real> &x) {
582 Real cnorm(0);
583 if (!isStateComputed_) {
584 int eflag(0);
585 const Real one(1);
586 const PartitionedVector<Real> &xp = partition(x);
587 // Set initial condition
588 uhist_[0]->set(*u0_);
589 if (useSketch_) { // Set initial guess for solve
590 stateSketch_->update();
591 uhist_[1]->set(*u0_);
592 }
593 // Run time stepper
594 size_type index;
595 for (size_type k = 1; k < Nt_; ++k) {
596 index = (useSketch_ ? 1 : k);
597 if (!useSketch_) { // Set initial guess for solve
598 uhist_[index]->set(*uhist_[index-1]);
599 }
600 // Solve state on current time interval
601 con_->update_uo(*uhist_[index-1], timeStamp_[k]);
602 con_->update_z(*xp.get(k), timeStamp_[k]);
603 con_->solve(*cprimal_, *uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
604 cnorm = std::max(cnorm,cprimal_->norm());
605 // Sketch state
606 if (useSketch_) {
607 eflag = stateSketch_->advance(one, *uhist_[1], static_cast<int>(k)-1, one);
608 throwError(eflag,"advance","solveState",574);
609 uhist_[0]->set(*uhist_[1]);
610 }
611 }
612 isStateComputed_ = true;
613 }
614 return cnorm;
615 }
616
617 Real updateSketch(const Vector<Real> &x, const Real tol) {
618 int eflag(0);
619 const PartitionedVector<Real> &xp = partition(x);
620 Real err(0), serr(0), cnorm(0); // cdot(0), dt(0);
621 bool flag = true;
622 while (flag) {
623 err = static_cast<Real>(0);
624 uhist_[0]->set(*u0_);
625 for (size_type k = 1; k < Nt_; ++k) {
626 eflag = stateSketch_->reconstruct(*uhist_[1],static_cast<int>(k)-1);
627 throwError(eflag,"reconstruct","updateSketch",592);
628 con_->update(*uhist_[0], *uhist_[1], *xp.get(k), timeStamp_[k]);
629 con_->value(*cprimal_, *uhist_[0], *uhist_[1], *xp.get(k), timeStamp_[k]);
630 /**** Linf norm for residual error ****/
631 cnorm = cprimal_->norm();
632 err = (cnorm > err ? cnorm : err);
633 if (err > tol) {
634 /**** L2 norm for residual error ****/
635 //cdot = cprimal_->dot(*cprimal_);
636 //dt = timeStamp_[k].t[1]-timeStamp_[k].t[0];
637 //err += cdot / dt;
638 //if (std::sqrt(err) > tol) {
639 break;
640 }
641 uhist_[0]->set(*uhist_[1]);
642 }
643 //err = std::sqrt(err);
644 if (stream_ != nullPtr) {
645 *stream_ << " *** State Rank: " << rankState_ << std::endl;
646 *stream_ << " *** Required Tolerance: " << tol << std::endl;
647 *stream_ << " *** Residual Norm: " << err << std::endl;
648 }
649 if (err > tol) {
650 //Real a(0.1838), b(3.1451); // Navier-Stokes
651 //Real a(2.6125), b(2.4841); // Semilinear
652 //rankState_ = std::max(rankState_+2,static_cast<size_t>(std::ceil((b-std::log(tol))/a)));
653 rankState_ *= updateFactor_; // Perhaps there is a better update strategy
655 stateSketch_->setRank(rankState_);
656 if (syncHessRank_) {
661 }
662 isStateComputed_ = false;
663 isAdjointComputed_ = false;
664 serr = solveState(x);
665 if (stream_ != nullPtr) {
666 *stream_ << " *** Maximum Solver Error: " << serr << std::endl;
667 }
668 }
669 else {
670 flag = false;
671 break;
672 }
673 }
674 return err;
675 }
676
677 /***************************************************************************/
678 /************ Methods to solve adjoint equation ****************************/
679 /***************************************************************************/
681 const Vector<Real> &uold, const Vector<Real> &unew,
682 const Vector<Real> &z, const TimeStamp<Real> &ts) {
683 obj_->gradient_un(*udual_, uold, unew, z, ts);
684 con_->applyInverseAdjointJacobian_un(l, *udual_, uold, unew, z, ts);
685 }
686
688 const Vector<Real> &uold, const Vector<Real> &unew,
689 const Vector<Real> &z, const TimeStamp<Real> &ts) {
690 const Real one(1);
691 obj_->gradient_uo(rhs, uold, unew, z, ts);
692 con_->applyAdjointJacobian_uo(*udual_, l, uold, unew, z, ts);
693 rhs.axpy(-one,*udual_);
694 }
695
697 const Vector<Real> &uold, const Vector<Real> &unew,
698 const Vector<Real> &z, const TimeStamp<Real> &ts) {
699 obj_->gradient_un(*udual_, uold, unew, z, ts);
700 rhs.plus(*udual_);
701 con_->applyInverseAdjointJacobian_un(l, rhs, uold, unew, z, ts);
702 }
703
704 void solveAdjoint(const Vector<Real> &x) {
705 int eflag(0);
706 if (!isAdjointComputed_) {
707 const PartitionedVector<Real> &xp = partition(x);
708 const Real one(1);
709 size_type uindex = (useSketch_ ? 1 : Nt_-1);
710 size_type lindex = (useSketch_ ? 0 : Nt_-1);
711 // Must first compute solve the state equation
712 solveState(x);
713 // Recover state from sketch
714 if (useSketch_) {
715 uhist_[1]->set(*uhist_[0]);
716 eflag = stateSketch_->reconstruct(*uhist_[0],static_cast<int>(Nt_)-3);
717 throwError(eflag,"reconstruct","solveAdjoint",672);
718 }
719 // Update dynamic constraint and objective
720 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
721 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
722 // Solve for terminal condition
724 *uhist_[uindex-1], *uhist_[uindex],
725 *xp.get(Nt_-1), timeStamp_[Nt_-1]);
726 if (useSketch_) {
727 eflag = adjointSketch_->advance(one,*lhist_[lindex],static_cast<int>(Nt_)-2,one);
728 throwError(eflag,"advance","solveAdjoint",683);
729 }
730 // Run reverse time stepper
731 for (size_type k = Nt_-2; k > 0; --k) {
732 // Compute k+1 component of rhs
733 computeAdjointRHS(*rhs_, *lhist_[lindex],
734 *uhist_[uindex-1], *uhist_[uindex],
735 *xp.get(k+1), timeStamp_[k+1]);
736 // Recover state from sketch
737 if (useSketch_) {
738 uhist_[1]->set(*uhist_[0]);
739 if (k==1) {
740 uhist_[0]->set(*u0_);
741 }
742 else {
743 eflag = stateSketch_->reconstruct(*uhist_[0],static_cast<int>(k)-2);
744 throwError(eflag,"reconstruct","solveAdjoint",699);
745 }
746 }
747 uindex = (useSketch_ ? 1 : k);
748 lindex = (useSketch_ ? 0 : k);
749 // Update dynamic constraint and objective
750 con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
751 obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
752 // Solve for adjoint on interval k
753 advanceAdjoint(*lhist_[lindex], *rhs_,
754 *uhist_[uindex-1], *uhist_[uindex],
755 *xp.get(k), timeStamp_[k]);
756 if (useSketch_) {
757 eflag = adjointSketch_->advance(one,*lhist_[lindex],static_cast<int>(k)-1,one);
758 throwError(eflag,"advance","solveAdjoint",713);
759 }
760 }
761 isAdjointComputed_ = true;
762 }
763 }
764
765 /***************************************************************************/
766 /************ Method for gradient computation ******************************/
767 /***************************************************************************/
769 const Vector<Real> &uold, const Vector<Real> &unew,
770 const Vector<Real> &z, const TimeStamp<Real> &ts) {
771 const Real one(1);
772 obj_->gradient_z(g, uold, unew, z, ts);
773 con_->applyAdjointJacobian_z(*zdual_, l, uold, unew, z, ts);
774 g.axpy(-one,*zdual_);
775 }
776
777 /***************************************************************************/
778 /************ Method to solve state sensitivity equation *******************/
779 /***************************************************************************/
781 const Vector<Real> &wold, const Vector<Real> &uold,
782 const Vector<Real> &unew, const Vector<Real> &z,
783 const TimeStamp<Real> &ts) {
784 const Real one(1);
785 con_->applyJacobian_z(*crhs_, v, uold, unew, z, ts);
786 con_->applyJacobian_uo(*cprimal_, wold, uold, unew, z, ts);
787 crhs_->axpy(-one, *cprimal_);
788 con_->applyInverseJacobian_un(wnew, *crhs_, uold, unew, z, ts);
789 }
790
791 /***************************************************************************/
792 /************ Methods to solve adjoint sensitivity equation ****************/
793 /***************************************************************************/
795 const Vector<Real> &v, const Vector<Real> &l,
796 const Vector<Real> &wold, const Vector<Real> &wnew,
797 const Vector<Real> &uold, const Vector<Real> &unew,
798 const Vector<Real> &z, const TimeStamp<Real> &ts) {
799 const Real one(1);
800 // Mixed derivative rhs term
801 con_->applyAdjointHessian_z_un(*rhs_, l, v, uold, unew, z, ts);
802 obj_->hessVec_un_z(*udual_, v, uold, unew, z, ts);
803 rhs_->axpy(-one, *udual_);
804 // State derivative rhs term
805 con_->applyAdjointHessian_un_un(*udual_, l, wnew, uold, unew, z, ts);
806 rhs_->axpy(-one, *udual_);
807 obj_->hessVec_un_un(*udual_, wnew, uold, unew, z, ts);
808 rhs_->plus(*udual_);
809 con_->applyAdjointHessian_uo_un(*udual_, l, wold, uold, unew, z, ts);
810 rhs_->axpy(-one, *udual_);
811 obj_->hessVec_un_uo(*udual_, wold, uold, unew, z, ts);
812 rhs_->plus(*udual_);
813 // Invert adjoint Jacobian
814 con_->applyInverseAdjointJacobian_un(p, *rhs_, uold, unew, z, ts);
815 }
816
818 const Vector<Real> &wold, const Vector<Real> &wnew,
819 const Vector<Real> &uold, const Vector<Real> &unew,
820 const Vector<Real> &z, const TimeStamp<Real> &ts,
821 const bool sumInto = false) {
822 const Real one(1);
823 // Compute new/old Hessian of Lagrangian
824 if (!sumInto) {
825 obj_->hessVec_un_uo(Hv, wold, uold, unew, z, ts);
826 }
827 else {
828 obj_->hessVec_un_uo(*udual_, wold, uold, unew, z, ts);
829 Hv.plus(*udual_);
830 }
831 con_->applyAdjointHessian_uo_un(*udual_, l, wold, uold, unew, z, ts);
832 Hv.axpy(-one,*udual_);
833 // Compute new/new Hessian of Lagrangian
834 obj_->hessVec_un_un(*udual_, wnew, uold, unew, z, ts);
835 Hv.plus(*udual_);
836 con_->applyAdjointHessian_un_un(*udual_, l, wnew, uold, unew, z, ts);
837 Hv.axpy(-one,*udual_);
838 }
839
841 const Vector<Real> &v, const Vector<Real> &l,
842 const Vector<Real> &uold, const Vector<Real> &unew,
843 const Vector<Real> &z, const TimeStamp<Real> &ts,
844 const bool sumInto = false) {
845 const Real one(1);
846 // Compute new/old Hessian of Lagrangian
847 if (!sumInto) {
848 con_->applyAdjointHessian_z_un(Hv, l, v, uold, unew, z, ts);
849 }
850 else {
851 con_->applyAdjointHessian_z_un(*udual_, l, v, uold, unew, z, ts);
852 Hv.plus(*udual_);
853 }
854 obj_->hessVec_un_z(*udual_, v, uold, unew, z, ts);
855 Hv.axpy(-one, *udual_);
856 }
857
859 const Vector<Real> &uold, const Vector<Real> &unew,
860 const Vector<Real> &z, const TimeStamp<Real> &ts,
861 const bool sumInto = false) {
862 const Real one(1);
863 if (!sumInto) {
864 con_->applyAdjointJacobian_uo(Hv, p, uold, unew, z, ts);
865 Hv.scale(-one);
866 }
867 else {
868 con_->applyAdjointJacobian_uo(*udual_, p, uold, unew, z, ts);
869 Hv.axpy(-one, *udual_);
870 }
871 }
872
874 const Vector<Real> &wold, const Vector<Real> &wnew,
875 const Vector<Real> &uold, const Vector<Real> &unew,
876 const Vector<Real> &z, const TimeStamp<Real> &ts,
877 const bool sumInto = false) {
878 const Real one(1);
879 // Compute old/new Hessian of Lagrangian
880 if (!sumInto) {
881 obj_->hessVec_uo_un(Hv, wnew, uold, unew, z, ts);
882 }
883 else {
884 obj_->hessVec_uo_un(*udual_, wnew, uold, unew, z, ts);
885 Hv.plus(*udual_);
886 }
887 con_->applyAdjointHessian_un_uo(*udual_, l, wnew, uold, unew, z, ts);
888 Hv.axpy(-one,*udual_);
889 // Compute old/old Hessian of Lagrangian
890 obj_->hessVec_uo_uo(*udual_, wold, uold, unew, z, ts);
891 Hv.plus(*udual_);
892 con_->applyAdjointHessian_uo_uo(*udual_, l, wold, uold, unew, z, ts);
893 Hv.axpy(-one,*udual_);
894 }
895
897 const Vector<Real> &v, const Vector<Real> &l,
898 const Vector<Real> &uold, const Vector<Real> &unew,
899 const Vector<Real> &z, const TimeStamp<Real> &ts,
900 const bool sumInto = false) {
901 const Real one(1);
902 // Compute new/old Hessian of Lagrangian
903 if (!sumInto) {
904 con_->applyAdjointHessian_z_uo(Hv, l, v, uold, unew, z, ts);
905 }
906 else {
907 con_->applyAdjointHessian_z_uo(*udual_, l, v, uold, unew, z, ts);
908 Hv.plus(*udual_);
909 }
910 obj_->hessVec_uo_z(*udual_, v, uold, unew, z, ts);
911 Hv.axpy(-one,*udual_);
912 }
913
915 const Vector<Real> &uold, const Vector<Real> &unew,
916 const Vector<Real> &z, const TimeStamp<Real> &ts) {
917 // Solve adjoint sensitivity on current time interval
918 con_->applyInverseAdjointJacobian_un(p, rhs, uold, unew, z, ts);
919 }
920
921 /***************************************************************************/
922 /************ Method for Hessian-times-a-vector computation ****************/
923 /***************************************************************************/
925 const Vector<Real> &v, const Vector<Real> &l,
926 const Vector<Real> &uold, const Vector<Real> &unew,
927 const Vector<Real> &z, const TimeStamp<Real> &ts) {
928 const Real one(1);
929 // Compute Hessian of Lagrangian
930 obj_->hessVec_z_z(Hv, v, uold, unew, z, ts);
931 con_->applyAdjointHessian_z_z(*zdual_, l, v, uold, unew, z, ts);
932 Hv.axpy(-one, *zdual_);
933 }
934
936 const Vector<Real> &wold, const Vector<Real> &wnew,
937 const Vector<Real> &uold, const Vector<Real> &unew,
938 const Vector<Real> &z, const TimeStamp<Real> &ts) {
939 const Real one(1);
940 // Compute Hessian of Lagrangian on previous time interval
941 obj_->hessVec_z_uo(*zdual_, wold, uold, unew, z, ts);
942 Hv.axpy(-one, *zdual_);
943 con_->applyAdjointHessian_uo_z(*zdual_, l, wold, uold, unew, z, ts);
944 Hv.plus(*zdual_);
945 // Compute Hessian of Lagrangian on current time interval
946 obj_->hessVec_z_un(*zdual_, wnew, uold, unew, z, ts);
947 Hv.axpy(-one, *zdual_);
948 con_->applyAdjointHessian_un_z(*zdual_, l, wnew, uold, unew, z, ts);
949 Hv.plus(*zdual_);
950 }
951
953 const Vector<Real> &uold, const Vector<Real> &unew,
954 const Vector<Real> &z, const TimeStamp<Real> &ts) {
955 con_->applyAdjointJacobian_z(*zdual_, p, uold, unew, z, ts);
956 Hv.plus(*zdual_);
957 }
958};
959
960} // namespace ROL
961
962#endif // ROL_REDUCEDDYNAMICOBJECTIVE_HPP
Defines the time-dependent constraint operator interface for simulation-based optimization.
Defines the time-dependent objective function interface for simulation-based optimization....
Provides the interface to evaluate objective functions.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Defines the linear algebra of vector space on a generic partitioned vector.
ROL::Ptr< const Vector< Real > > get(size_type i) const
static Ptr< PartitionedVector > create(std::initializer_list< Vp > vs)
Defines the reduced time-dependent objective function interface for simulation-based optimization.
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Real updateSketch(const Vector< Real > &x, const Real tol)
const PartitionedVector< Real > & partition(const Vector< Real > &x) const
Ptr< Vector< Real > > makeDynamicVector(const Vector< Real > &x) const
void addMixedHessLag(Vector< Real > &Hv, const Vector< Real > &l, const Vector< Real > &wold, const Vector< Real > &wnew, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
void advanceStateSens(Vector< Real > &wnew, const Vector< Real > &v, const Vector< Real > &wold, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
void computeNewMixedHessLag(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts, const bool sumInto=false)
void computeNewStateJacobian(Vector< Real > &Hv, const Vector< Real > &p, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts, const bool sumInto=false)
typename std::vector< Real >::size_type size_type
void addAdjointSens(Vector< Real > &Hv, const Vector< Real > &p, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
void setTerminalConditionHess(Vector< Real > &p, const Vector< Real > &v, const Vector< Real > &l, const Vector< Real > &wold, const Vector< Real > &wnew, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
void computeOldMixedHessLag(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts, const bool sumInto=false)
std::vector< Ptr< Vector< Real > > > uhist_
std::vector< Ptr< Vector< Real > > > phist_
Real solveState(const Vector< Real > &x)
void computeNewStateHessLag(Vector< Real > &Hv, const Vector< Real > &l, const Vector< Real > &wold, const Vector< Real > &wnew, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts, const bool sumInto=false)
void setTerminalCondition(Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
std::vector< Ptr< Vector< Real > > > whist_
const std::vector< TimeStamp< Real > > timeStamp_
void solveAdjoint(const Vector< Real > &x)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
void updateGradient(Vector< Real > &g, const Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
Real value(const Vector< Real > &x, Real &tol)
Compute value.
void advanceAdjointSens(Vector< Real > &p, Vector< Real > &rhs, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
void computeOldStateHessLag(Vector< Real > &Hv, const Vector< Real > &l, const Vector< Real > &wold, const Vector< Real > &wnew, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts, const bool sumInto=false)
const Ptr< DynamicConstraint< Real > > con_
std::vector< Ptr< Vector< Real > > > lhist_
void advanceAdjoint(Vector< Real > &l, Vector< Real > &rhs, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
PartitionedVector< Real > & partition(Vector< Real > &x) const
const Ptr< DynamicObjective< Real > > obj_
void computeControlHessLag(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
ReducedDynamicObjective(const Ptr< DynamicObjective< Real > > &obj, const Ptr< DynamicConstraint< Real > > &con, const Ptr< Vector< Real > > &u0, const Ptr< Vector< Real > > &zvec, const Ptr< Vector< Real > > &cvec, const std::vector< TimeStamp< Real > > &timeStamp, ROL::ParameterList &pl, const Ptr< std::ostream > &stream=nullPtr)
void computeAdjointRHS(Vector< Real > &rhs, const Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
void throwError(const int err, const std::string &sfunc, const std::string &func, const int line) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void scale(const Real alpha)=0
Compute where .
virtual void print(std::ostream &outStream) const
Definition: ROL_Vector.hpp:261
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
Contains local time step information.