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"
50 #include "ROL_DynamicObjective.hpp"
52 
53 
76 namespace ROL {
77 
78 template<typename Real>
79 class ReducedDynamicObjective : public Objective<Real> {
81 private:
82  // Problem data.
83  const Ptr<DynamicObjective<Real>> obj_;
84  const Ptr<DynamicConstraint<Real>> con_;
85  const Ptr<Vector<Real>> u0_;
86  const std::vector<TimeStamp<Real>> timeStamp_;
87  const size_type Nt_;
88  // General sketch information.
89  const bool useSketch_;
90  // State sketch information.
92  Ptr<Sketch<Real>> stateSketch_;
93  // Adjoint sketch information.
94  const int rankAdjoint_;
95  Ptr<Sketch<Real>> adjointSketch_;
96  // State sensitivity sketch information.
97  const int rankStateSens_;
98  Ptr<Sketch<Real>> stateSensSketch_;
99  // Inexactness information.
100  const bool useInexact_;
101  const Real updateFactor_;
102  const int maxRank_;
103  // Vector storage for intermediate computations.
104  std::vector<Ptr<Vector<Real>>> uhist_; // State history.
105  std::vector<Ptr<Vector<Real>>> lhist_; // Adjoint history.
106  std::vector<Ptr<Vector<Real>>> whist_; // State sensitivity history.
107  std::vector<Ptr<Vector<Real>>> phist_; // Adjoint sensitivity history.
108  Ptr<Vector<Real>> cprimal_; // Primal constraint vector.
109  Ptr<Vector<Real>> crhs_; // Primal constraint vector.
110  Ptr<Vector<Real>> udual_; // Dual state vector.
111  Ptr<Vector<Real>> rhs_; // Dual state vector.
112  Ptr<Vector<Real>> zdual_; // Dual control vector.
113  Real val_; // Objective function value.
114  // Flags.
115  bool isValueComputed_; // Whether objective has been evaluated.
116  bool isStateComputed_; // Whether state has been solved.
117  bool isAdjointComputed_; // Whether adjoint has been solved.
118  bool useHessian_; // Whether to use Hessian or FD.
119  bool useSymHess_; // Whether to use symmetric Hessian approximation.
120 
122  return static_cast<PartitionedVector<Real>&>(x);
123  }
124 
125  const PartitionedVector<Real>& partition ( const Vector<Real>& x ) const {
126  return static_cast<const PartitionedVector<Real>&>(x);
127  }
128 
129 public:
131  const Ptr<DynamicConstraint<Real>> &con,
132  const Ptr<Vector<Real>> &u0,
133  const Ptr<Vector<Real>> &zvec,
134  const Ptr<Vector<Real>> &cvec,
135  const std::vector<TimeStamp<Real>> &timeStamp,
136  ROL::ParameterList &pl)
137  : obj_ ( obj ), // Dynamic objective function.
138  con_ ( con ), // Dynamic constraint function.
139  u0_ ( u0 ), // Initial condition.
140  timeStamp_ ( timeStamp ), // Vector of time stamps.
141  Nt_ ( timeStamp.size() ), // Number of time intervals.
142  useSketch_ ( pl.get("Use Sketching", true) ), // Use state sketch if true.
143  rankState_ ( pl.get("State Rank", 10) ), // Rank of state sketch.
144  stateSketch_ ( nullPtr ), // State sketch object.
145  rankAdjoint_ ( pl.get("Adjoint Rank", 10) ), // Rank of adjoint sketch.
146  adjointSketch_ ( nullPtr ), // Adjoint sketch object.
147  rankStateSens_ ( pl.get("State Sensitivity Rank", 10) ), // Rank of state sensitivity sketch.
148  stateSensSketch_ ( nullPtr ), // State sensitivity sketch object.
149  useInexact_ ( pl.get("Adaptive Rank", false) ), // Update rank adaptively.
150  updateFactor_ ( pl.get("Rank Update Factor", 2.0) ), // Rank update factor.
151  maxRank_ ( pl.get("Maximum Rank", Nt_) ), // Maximum rank.
152  isValueComputed_ ( false ), // Flag indicating whether value has been computed.
153  isStateComputed_ ( false ), // Flag indicating whether state has been computed.
154  isAdjointComputed_ ( false ), // Flag indicating whether adjoint has been computed.
155  useHessian_ ( pl.get("Use Hessian", true) ), // Flag indicating whether to use the Hessian.
156  useSymHess_ ( pl.get("Use Only Sketched Sensitivity", true) ) {// Flag indicating whether to use symmetric sketched Hessian.
157  uhist_.clear(); lhist_.clear(); whist_.clear(); phist_.clear();
158  if (useSketch_) { // Only maintain a sketch of the state time history
159  stateSketch_ = makePtr<Sketch<Real>>(*u0_,static_cast<int>(Nt_)-1,rankState_);
160  uhist_.push_back(u0_->clone());
161  uhist_.push_back(u0_->clone());
162  lhist_.push_back(cvec->dual().clone());
163  adjointSketch_ = makePtr<Sketch<Real>>(*u0_,static_cast<int>(Nt_)-1,rankAdjoint_);
164  if (useHessian_) {
165  stateSensSketch_ = makePtr<Sketch<Real>>(*u0_,static_cast<int>(Nt_)-1,rankStateSens_);
166  whist_.push_back(u0_->clone());
167  whist_.push_back(u0_->clone());
168  phist_.push_back(cvec->dual().clone());
169  }
170  }
171  else { // Store entire state time history
172  for (size_type i = 0; i < Nt_; ++i) {
173  uhist_.push_back(u0_->clone());
174  lhist_.push_back(cvec->dual().clone());
175  if (useHessian_) {
176  whist_.push_back(u0_->clone());
177  phist_.push_back(cvec->dual().clone());
178  }
179  }
180  }
181  cprimal_ = cvec->clone();
182  crhs_ = cvec->clone();
183  udual_ = u0_->dual().clone();
184  rhs_ = u0_->dual().clone();
185  zdual_ = zvec->dual().clone();
186  }
187 
188  Ptr<Vector<Real>> makeDynamicVector(const Vector<Real> &x) const {
190  }
191 
192  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
193  if (useSketch_) {
194  stateSketch_->update();
195  if (flag == true) {
196  adjointSketch_->update();
197  if (useHessian_) {
198  stateSensSketch_->update();
199  }
200  }
201  }
202  for (size_type i = 0; i < uhist_.size(); ++i) {
203  uhist_[i]->zero();
204  }
205  val_ = static_cast<Real>(0);
206  isValueComputed_ = false;
207  isStateComputed_ = false;
208  if (flag == true) {
209  isAdjointComputed_ = false;
210  }
211  }
212 
213  Real value( const Vector<Real> &x, Real &tol ) {
214  if (!isValueComputed_) {
215  const Real one(1);
216  const PartitionedVector<Real> &xp = partition(x);
217  // Set initial condition
218  uhist_[0]->set(*u0_);
219  if (useSketch_) {
220  stateSketch_->update();
221  //stateSketch_->advance(one,*uhist_[0],0,one);
222  }
223  // Run time stepper
224  Real valk(0);
225  size_type index;
226  for (size_type k = 1; k < Nt_; ++k) {
227  index = (useSketch_ ? 1 : k);
228  // Update dynamic constraint and objective
229  con_->update(*uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
230  obj_->update(*uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
231  // Solve state on current time interval
232  con_->solve(*cprimal_, *uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
233  // Compute objective function value on current time interval
234  valk = obj_->value(*uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
235  // Update total objective function value
236  val_ += valk;
237  // Sketch state
238  if (useSketch_) {
239  stateSketch_->advance(one,*uhist_[1],static_cast<int>(k)-1,one);
240  uhist_[0]->set(*uhist_[1]);
241  }
242  }
243  isValueComputed_ = true;
244  isStateComputed_ = true;
245  }
246  return val_;
247  }
248 
249  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
251  const PartitionedVector<Real> &xp = partition(x);
252  const Real one(1);
253  size_type uindex = (useSketch_ ? 1 : Nt_-1);
254  size_type lindex = (useSketch_ ? 0 : Nt_-1);
255  // Must first compute the value
256  solveState(x);
257  if (useSketch_ && useInexact_) {
258  tol = updateSketch(x,tol);
259  }
260  // Recover state from sketch
261  if (useSketch_) {
262  uhist_[1]->set(*uhist_[0]);
263  stateSketch_->reconstruct(*uhist_[0],static_cast<int>(Nt_)-3);
264  if (isAdjointComputed_) {
265  adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(Nt_)-2);
266  }
267  }
268  // Update dynamic constraint and objective
269  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
270  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
271  // Solve for terminal condition
272  if (!isAdjointComputed_) {
273  setTerminalCondition(*lhist_[lindex],
274  *uhist_[uindex-1], *uhist_[uindex],
275  *xp.get(Nt_-1), timeStamp_[Nt_-1]);
276  if (useSketch_) {
277  adjointSketch_->advance(one,*lhist_[0],static_cast<int>(Nt_)-2,one);
278  }
279  }
280  // Update gradient on terminal interval
281  updateGradient(*gp.get(Nt_-1), *lhist_[lindex],
282  *uhist_[uindex-1], *uhist_[uindex],
283  *xp.get(Nt_-1), timeStamp_[Nt_-1]);
284  // Run reverse time stepper
285  for (size_type k = Nt_-2; k > 0; --k) {
286  if (!isAdjointComputed_) {
287  // Compute k+1 component of rhs
288  computeAdjointRHS(*rhs_, *lhist_[lindex],
289  *uhist_[uindex-1], *uhist_[uindex],
290  *xp.get(k+1), timeStamp_[k+1]);
291  }
292  uindex = (useSketch_ ? 1 : k);
293  lindex = (useSketch_ ? 0 : k);
294  // Recover state from sketch
295  if (useSketch_) {
296  uhist_[1]->set(*uhist_[0]);
297  if (k==1) {
298  uhist_[0]->set(*u0_);
299  }
300  else {
301  stateSketch_->reconstruct(*uhist_[0],static_cast<int>(k)-2);
302  }
303  if (isAdjointComputed_) {
304  adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(k)-1);
305  }
306  }
307  // Update dynamic constraint and objective
308  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
309  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
310  // Solve for adjoint on interval k
311  if (!isAdjointComputed_) {
312  advanceAdjoint(*lhist_[lindex], *rhs_,
313  *uhist_[uindex-1], *uhist_[uindex],
314  *xp.get(k), timeStamp_[k]);
315  if (useSketch_) {
316  adjointSketch_->advance(one,*lhist_[0],static_cast<int>(k)-1,one);
317  }
318  }
319  // Update gradient on interval k
320  updateGradient(*gp.get(k), *lhist_[lindex],
321  *uhist_[uindex-1], *uhist_[uindex],
322  *xp.get(k), timeStamp_[k]);
323  }
324  isAdjointComputed_ = true;
325  }
326 
327  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
328  if (useHessian_) {
329  // Must first solve the state and adjoint equations
330  solveState(x);
331  solveAdjoint(x);
332  // Now compute Hessian
333  const Real one(1);
334  const PartitionedVector<Real> &xp = partition(x);
335  const PartitionedVector<Real> &vp = partition(v);
337  // Compute state sensitivity
338  whist_[0]->zero();
339  if (useSketch_) {
340  stateSensSketch_->update();
341  //stateSensSketch_->advance(one,*whist_[0],0,one);
342  uhist_[0]->set(*u0_);
343  }
344  size_type uindex, lindex;
345  for (size_type k = 1; k < Nt_; ++k) {
346  uindex = (useSketch_ ? 1 : k);
347  lindex = (useSketch_ ? 0 : k);
348  // Reconstruct sketched state
349  if (useSketch_) {
350  stateSketch_->reconstruct(*uhist_[1],static_cast<int>(k)-1);
351  adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(k)-1);
352  }
353  // Update dynamic constraint and objective
354  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
355  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
356  // Advance state sensitivity on current time interval
357  advanceStateSens(*whist_[uindex],
358  *vp.get(k), *whist_[uindex-1],
359  *uhist_[uindex-1], *uhist_[uindex],
360  *xp.get(k), timeStamp_[k]);
361  // Set control only Hessian
362  computeControlHessLag(*hvp.get(k),
363  *vp.get(k), *lhist_[lindex],
364  *uhist_[uindex-1], *uhist_[uindex],
365  *xp.get(k), timeStamp_[k]);
366  if (!useSymHess_) {
367  // Add mixed derivative Hessian
368  addMixedHessLag(*hvp.get(k), *lhist_[lindex],
369  *whist_[uindex-1], *whist_[uindex],
370  *uhist_[uindex-1], *uhist_[uindex],
371  *xp.get(k), timeStamp_[k]);
372  }
373  // Sketch state sensitivity
374  if (useSketch_) {
375  stateSensSketch_->advance(one,*whist_[1],static_cast<int>(k)-1,one);
376  whist_[0]->set(*whist_[1]);
377  uhist_[0]->set(*uhist_[1]);
378  }
379  }
380 
381  // Compute adjoint sensitivity
382  uindex = (useSketch_ ? 1 : Nt_-1);
383  lindex = (useSketch_ ? 0 : Nt_-1);
384  if (useSketch_) {
385  // Recover terminal state
386  uhist_[1]->set(*uhist_[0]);
387  stateSketch_->reconstruct(*uhist_[0],static_cast<int>(Nt_)-3);
388  // Recover terminal adjoint
389  adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(Nt_)-2);
390  // Recover terminal state sensitivity
391  whist_[1]->set(*whist_[0]);
392  stateSensSketch_->reconstruct(*whist_[0],static_cast<int>(Nt_)-3);
393  }
394  // Update dynamic constraint and objective
395  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
396  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
397  // Solve for terminal condition
399  *vp.get(Nt_-1), *lhist_[lindex],
400  *whist_[uindex-1], *whist_[uindex],
401  *uhist_[uindex-1], *uhist_[uindex],
402  *xp.get(Nt_-1), timeStamp_[Nt_-1]);
403  if (useSymHess_) {
404  // Add mixed derivative Hessian
405  addMixedHessLag(*hvp.get(Nt_-1), *lhist_[lindex],
406  *whist_[uindex-1], *whist_[uindex],
407  *uhist_[uindex-1], *uhist_[uindex],
408  *xp.get(Nt_-1), timeStamp_[Nt_-1]);
409  }
410  // Add adjoint sensitivity to Hessian
411  addAdjointSens(*hvp.get(Nt_-1), *phist_[lindex],
412  *uhist_[uindex-1], *uhist_[uindex],
413  *xp.get(Nt_-1), timeStamp_[Nt_-1]);
414  // Run reverse time stepper
415  for (size_type k = Nt_-2; k > 0; --k) {
416  // Compute new components of rhs
418  *whist_[uindex-1], *whist_[uindex],
419  *uhist_[uindex-1], *uhist_[uindex],
420  *xp.get(k+1), timeStamp_[k+1],
421  false);
423  *vp.get(k+1), *lhist_[lindex],
424  *uhist_[uindex-1], *uhist_[uindex],
425  *xp.get(k+1), timeStamp_[k+1],
426  true);
428  *uhist_[uindex-1], *uhist_[uindex],
429  *xp.get(k+1), timeStamp_[k+1],
430  true);
431  // Recover state, adjoint and state sensitivity from sketch
432  if (useSketch_) {
433  uhist_[1]->set(*uhist_[0]);
434  whist_[1]->set(*whist_[0]);
435  if (k==1) {
436  uhist_[0]->set(*u0_);
437  whist_[0]->zero();
438  }
439  else {
440  stateSketch_->reconstruct(*uhist_[0],static_cast<int>(k)-2);
441  stateSensSketch_->reconstruct(*whist_[0],static_cast<int>(k)-2);
442  }
443  adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(k)-1);
444  }
445  uindex = (useSketch_ ? 1 : k);
446  lindex = (useSketch_ ? 0 : k);
447  // Update dynamic constraint and objective
448  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
449  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
450  // Compute old components of rhs
452  *whist_[uindex-1], *whist_[uindex],
453  *uhist_[uindex-1], *uhist_[uindex],
454  *xp.get(k), timeStamp_[k],
455  true);
457  *vp.get(k), *lhist_[lindex],
458  *uhist_[uindex-1], *uhist_[uindex],
459  *xp.get(k), timeStamp_[k],
460  true);
461  // Solve for adjoint on interval k
462  advanceAdjointSens(*phist_[lindex], *rhs_,
463  *uhist_[uindex-1], *uhist_[uindex],
464  *xp.get(k), timeStamp_[k]);
465  if (useSymHess_) {
466  // Add mixed derivative Hessian
467  addMixedHessLag(*hvp.get(k), *lhist_[lindex],
468  *whist_[uindex-1], *whist_[uindex],
469  *uhist_[uindex-1], *uhist_[uindex],
470  *xp.get(k), timeStamp_[k]);
471  }
472  // Add adjoint sensitivity to Hessian
473  addAdjointSens(*hvp.get(k), *phist_[lindex],
474  *uhist_[uindex-1], *uhist_[uindex],
475  *xp.get(k), timeStamp_[k]);
476  }
477  }
478  else {
479  Objective<Real>::hessVec(hv,v,x,tol);
480  }
481  }
482 
483 private:
484  /***************************************************************************/
485  /************ Method to solve state equation *******************************/
486  /***************************************************************************/
487  void solveState(const Vector<Real> &x) {
488  if (!isStateComputed_) {
489  const Real one(1);
490  const PartitionedVector<Real> &xp = partition(x);
491  // Set initial condition
492  uhist_[0]->set(*u0_);
493  if (useSketch_) {
494  //stateSketch_->advance(one,*uhist_[0],0,one);
495  }
496  // Run time stepper
497  size_type index;
498  for (size_type k = 1; k < Nt_; ++k) {
499  index = (useSketch_ ? 1 : k);
500  // Update dynamic constraint and objective
501  con_->update(*uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
502  obj_->update(*uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
503  // Solve state on current time interval
504  con_->solve(*cprimal_, *uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
505  // Sketch state
506  if (useSketch_) {
507  stateSketch_->advance(one,*uhist_[index],static_cast<int>(k)-1,one);
508  uhist_[index-1]->set(*uhist_[index]);
509  }
510  }
511  isStateComputed_ = true;
512  }
513  }
514 
515  Real updateSketch(const Vector<Real> &x, const Real tol) {
516  const PartitionedVector<Real> &xp = partition(x);
517  Real err(0), cnorm(0);
518  bool flag = true;
519  while (flag) {
520  err = static_cast<Real>(0);
521  uhist_[0]->set(*u0_);
522  for (size_type k = 1; k < Nt_; ++k) {
523  stateSketch_->reconstruct(*uhist_[1],static_cast<int>(k));
524  con_->value(*cprimal_, *uhist_[0], *uhist_[1], *xp.get(k), timeStamp_[k]);
525  cnorm = cprimal_->norm();
526  err = (cnorm > err ? cnorm : err); // Use Linf norm. Could change to use, e.g., L2.
527  if (err > tol) {
528  break;
529  }
530  }
531  if (err > tol) {
532  rankState_ *= updateFactor_; // Perhaps there is a better update strategy
534  stateSketch_->setRank(rankState_);
535  solveState(x);
536  }
537  else {
538  flag = false;
539  break;
540  }
541  }
542  return err;
543  }
544 
545  /***************************************************************************/
546  /************ Methods to solve adjoint equation ****************************/
547  /***************************************************************************/
549  const Vector<Real> &uold, const Vector<Real> &unew,
550  const Vector<Real> &z, const TimeStamp<Real> &ts) {
551  obj_->gradient_un(*udual_, uold, unew, z, ts);
552  con_->applyInverseAdjointJacobian_un(l, *udual_, uold, unew, z, ts);
553  }
554 
556  const Vector<Real> &uold, const Vector<Real> &unew,
557  const Vector<Real> &z, const TimeStamp<Real> &ts) {
558  const Real one(1);
559  obj_->gradient_uo(rhs, uold, unew, z, ts);
560  con_->applyAdjointJacobian_uo(*udual_, l, uold, unew, z, ts);
561  rhs.axpy(-one,*udual_);
562  }
563 
565  const Vector<Real> &uold, const Vector<Real> &unew,
566  const Vector<Real> &z, const TimeStamp<Real> &ts) {
567  obj_->gradient_un(*udual_, uold, unew, z, ts);
568  rhs.plus(*udual_);
569  con_->applyInverseAdjointJacobian_un(l, rhs, uold, unew, z, ts);
570  }
571 
572  void solveAdjoint(const Vector<Real> &x) {
573  if (!isAdjointComputed_) {
574  const PartitionedVector<Real> &xp = partition(x);
575  const Real one(1);
576  size_type uindex = (useSketch_ ? 1 : Nt_-1);
577  size_type lindex = (useSketch_ ? 0 : Nt_-1);
578  // Must first compute solve the state equation
579  solveState(x);
580  // Recover state from sketch
581  if (useSketch_) {
582  uhist_[1]->set(*uhist_[0]);
583  stateSketch_->reconstruct(*uhist_[0],static_cast<int>(Nt_)-3);
584  }
585  // Update dynamic constraint and objective
586  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
587  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
588  // Solve for terminal condition
589  setTerminalCondition(*lhist_[lindex],
590  *uhist_[uindex-1], *uhist_[uindex],
591  *xp.get(Nt_-1), timeStamp_[Nt_-1]);
592  if (useSketch_) {
593  adjointSketch_->advance(one,*lhist_[lindex],static_cast<int>(Nt_)-2,one);
594  }
595  // Run reverse time stepper
596  for (size_type k = Nt_-2; k > 0; --k) {
597  // Compute k+1 component of rhs
598  computeAdjointRHS(*rhs_, *lhist_[lindex],
599  *uhist_[uindex-1], *uhist_[uindex],
600  *xp.get(k+1), timeStamp_[k+1]);
601  // Recover state from sketch
602  if (useSketch_) {
603  uhist_[1]->set(*uhist_[0]);
604  if (k==1) {
605  uhist_[0]->set(*u0_);
606  }
607  else {
608  stateSketch_->reconstruct(*uhist_[0],static_cast<int>(k)-2);
609  }
610  }
611  uindex = (useSketch_ ? 1 : k);
612  lindex = (useSketch_ ? 0 : k);
613  // Update dynamic constraint and objective
614  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
615  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
616  // Solve for adjoint on interval k
617  advanceAdjoint(*lhist_[lindex], *rhs_,
618  *uhist_[uindex-1], *uhist_[uindex],
619  *xp.get(k), timeStamp_[k]);
620  if (useSketch_) {
621  adjointSketch_->advance(one,*lhist_[lindex],static_cast<int>(k)-1,one);
622  }
623  }
624  isAdjointComputed_ = true;
625  }
626  }
627 
628  /***************************************************************************/
629  /************ Method for gradient computation ******************************/
630  /***************************************************************************/
632  const Vector<Real> &uold, const Vector<Real> &unew,
633  const Vector<Real> &z, const TimeStamp<Real> &ts) {
634  const Real one(1);
635  obj_->gradient_z(g, uold, unew, z, ts);
636  con_->applyAdjointJacobian_z(*zdual_, l, uold, unew, z, ts);
637  g.axpy(-one,*zdual_);
638  }
639 
640  /***************************************************************************/
641  /************ Method to solve state sensitivity equation *******************/
642  /***************************************************************************/
644  const Vector<Real> &wold, const Vector<Real> &uold,
645  const Vector<Real> &unew, const Vector<Real> &z,
646  const TimeStamp<Real> &ts) {
647  const Real one(1);
648  con_->applyJacobian_z(*crhs_, v, uold, unew, z, ts);
649  con_->applyJacobian_uo(*cprimal_, wold, uold, unew, z, ts);
650  crhs_->axpy(-one, *cprimal_);
651  con_->applyInverseJacobian_un(wnew, *crhs_, uold, unew, z, ts);
652  }
653 
654  /***************************************************************************/
655  /************ Methods to solve adjoint sensitivity equation ****************/
656  /***************************************************************************/
658  const Vector<Real> &v, const Vector<Real> &l,
659  const Vector<Real> &wold, const Vector<Real> &wnew,
660  const Vector<Real> &uold, const Vector<Real> &unew,
661  const Vector<Real> &z, const TimeStamp<Real> &ts) {
662  const Real one(1);
663  // Mixed derivative rhs term
664  con_->applyAdjointHessian_z_un(*rhs_, l, v, uold, unew, z, ts);
665  obj_->hessVec_un_z(*udual_, v, uold, unew, z, ts);
666  rhs_->axpy(-one, *udual_);
667  // State derivative rhs term
668  con_->applyAdjointHessian_un_un(*udual_, l, wnew, uold, unew, z, ts);
669  rhs_->axpy(-one, *udual_);
670  obj_->hessVec_un_un(*udual_, wnew, uold, unew, z, ts);
671  rhs_->plus(*udual_);
672  con_->applyAdjointHessian_uo_un(*udual_, l, wold, uold, unew, z, ts);
673  rhs_->axpy(-one, *udual_);
674  obj_->hessVec_un_uo(*udual_, wold, uold, unew, z, ts);
675  rhs_->plus(*udual_);
676  // Invert adjoint Jacobian
677  con_->applyInverseAdjointJacobian_un(p, *rhs_, uold, unew, z, ts);
678  }
679 
681  const Vector<Real> &wold, const Vector<Real> &wnew,
682  const Vector<Real> &uold, const Vector<Real> &unew,
683  const Vector<Real> &z, const TimeStamp<Real> &ts,
684  const bool sumInto = false) {
685  const Real one(1);
686  // Compute new/old Hessian of Lagrangian
687  if (!sumInto) {
688  obj_->hessVec_un_uo(Hv, wold, uold, unew, z, ts);
689  }
690  else {
691  obj_->hessVec_un_uo(*udual_, wold, uold, unew, z, ts);
692  Hv.plus(*udual_);
693  }
694  con_->applyAdjointHessian_uo_un(*udual_, l, wold, uold, unew, z, ts);
695  Hv.axpy(-one,*udual_);
696  // Compute new/new Hessian of Lagrangian
697  obj_->hessVec_un_un(*udual_, wnew, uold, unew, z, ts);
698  Hv.plus(*udual_);
699  con_->applyAdjointHessian_un_un(*udual_, l, wnew, uold, unew, z, ts);
700  Hv.axpy(-one,*udual_);
701  }
702 
704  const Vector<Real> &v, const Vector<Real> &l,
705  const Vector<Real> &uold, const Vector<Real> &unew,
706  const Vector<Real> &z, const TimeStamp<Real> &ts,
707  const bool sumInto = false) {
708  const Real one(1);
709  // Compute new/old Hessian of Lagrangian
710  if (!sumInto) {
711  con_->applyAdjointHessian_z_un(Hv, l, v, uold, unew, z, ts);
712  }
713  else {
714  con_->applyAdjointHessian_z_un(*udual_, l, v, uold, unew, z, ts);
715  Hv.plus(*udual_);
716  }
717  obj_->hessVec_un_z(*udual_, v, uold, unew, z, ts);
718  Hv.axpy(-one, *udual_);
719  }
720 
722  const Vector<Real> &uold, const Vector<Real> &unew,
723  const Vector<Real> &z, const TimeStamp<Real> &ts,
724  const bool sumInto = false) {
725  const Real one(1);
726  if (!sumInto) {
727  con_->applyAdjointJacobian_uo(Hv, p, uold, unew, z, ts);
728  Hv.scale(-one);
729  }
730  else {
731  con_->applyAdjointJacobian_uo(*udual_, p, uold, unew, z, ts);
732  Hv.axpy(-one, *udual_);
733  }
734  }
735 
737  const Vector<Real> &wold, const Vector<Real> &wnew,
738  const Vector<Real> &uold, const Vector<Real> &unew,
739  const Vector<Real> &z, const TimeStamp<Real> &ts,
740  const bool sumInto = false) {
741  const Real one(1);
742  // Compute old/new Hessian of Lagrangian
743  if (!sumInto) {
744  obj_->hessVec_uo_un(Hv, wnew, uold, unew, z, ts);
745  }
746  else {
747  obj_->hessVec_uo_un(*udual_, wnew, uold, unew, z, ts);
748  Hv.plus(*udual_);
749  }
750  con_->applyAdjointHessian_un_uo(*udual_, l, wnew, uold, unew, z, ts);
751  Hv.axpy(-one,*udual_);
752  // Compute old/old Hessian of Lagrangian
753  obj_->hessVec_uo_uo(*udual_, wold, uold, unew, z, ts);
754  Hv.plus(*udual_);
755  con_->applyAdjointHessian_uo_uo(*udual_, l, wold, uold, unew, z, ts);
756  Hv.axpy(-one,*udual_);
757  }
758 
760  const Vector<Real> &v, const Vector<Real> &l,
761  const Vector<Real> &uold, const Vector<Real> &unew,
762  const Vector<Real> &z, const TimeStamp<Real> &ts,
763  const bool sumInto = false) {
764  const Real one(1);
765  // Compute new/old Hessian of Lagrangian
766  if (!sumInto) {
767  con_->applyAdjointHessian_z_uo(Hv, l, v, uold, unew, z, ts);
768  }
769  else {
770  con_->applyAdjointHessian_z_uo(*udual_, l, v, uold, unew, z, ts);
771  Hv.plus(*udual_);
772  }
773  obj_->hessVec_uo_z(*udual_, v, uold, unew, z, ts);
774  Hv.axpy(-one,*udual_);
775  }
776 
778  const Vector<Real> &uold, const Vector<Real> &unew,
779  const Vector<Real> &z, const TimeStamp<Real> &ts) {
780  // Solve adjoint sensitivity on current time interval
781  con_->applyInverseAdjointJacobian_un(p, rhs, uold, unew, z, ts);
782  }
783 
784  /***************************************************************************/
785  /************ Method for Hessian-times-a-vector computation ****************/
786  /***************************************************************************/
788  const Vector<Real> &v, const Vector<Real> &l,
789  const Vector<Real> &uold, const Vector<Real> &unew,
790  const Vector<Real> &z, const TimeStamp<Real> &ts) {
791  const Real one(1);
792  // Compute Hessian of Lagrangian
793  obj_->hessVec_z_z(Hv, v, uold, unew, z, ts);
794  con_->applyAdjointHessian_z_z(*zdual_, l, v, uold, unew, z, ts);
795  Hv.axpy(-one, *zdual_);
796  }
797 
799  const Vector<Real> &wold, const Vector<Real> &wnew,
800  const Vector<Real> &uold, const Vector<Real> &unew,
801  const Vector<Real> &z, const TimeStamp<Real> &ts) {
802  const Real one(1);
803  // Compute Hessian of Lagrangian on previous time interval
804  obj_->hessVec_z_uo(*zdual_, wold, uold, unew, z, ts);
805  Hv.axpy(-one, *zdual_);
806  con_->applyAdjointHessian_uo_z(*zdual_, l, wold, uold, unew, z, ts);
807  Hv.plus(*zdual_);
808  // Compute Hessian of Lagrangian on current time interval
809  obj_->hessVec_z_un(*zdual_, wnew, uold, unew, z, ts);
810  Hv.axpy(-one, *zdual_);
811  con_->applyAdjointHessian_un_z(*zdual_, l, wnew, uold, unew, z, ts);
812  Hv.plus(*zdual_);
813  }
814 
816  const Vector<Real> &uold, const Vector<Real> &unew,
817  const Vector<Real> &z, const TimeStamp<Real> &ts) {
818  con_->applyAdjointJacobian_z(*zdual_, p, uold, unew, z, ts);
819  Hv.plus(*zdual_);
820  }
821 };
822 
823 } // namespace ROL
824 
825 #endif // ROL_REDUCEDDYNAMICOBJECTIVE_HPP
void advanceAdjointSens(Vector< Real > &p, Vector< Real > &rhs, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
Provides the interface to evaluate objective functions.
static Ptr< PartitionedVector > create(std::initializer_list< Vp > vs)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
typename PV< Real >::size_type size_type
virtual void scale(const Real alpha)=0
Compute where .
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)
Defines the reduced time-dependent objective function interface for simulation-based optimization...
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
Defines the time-dependent constraint operator interface for simulation-based optimization.
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)
Defines the linear algebra of vector space on a generic partitioned vector.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
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)
const std::vector< TimeStamp< Real > > timeStamp_
Contains local time step information.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
Defines the time-dependent objective function interface for simulation-based optimization. Computes time-local contributions of value, gradient, Hessian-vector product etc to a larger composite objective defined over the simulation time. In contrast to other objective classes Objective_TimeSimOpt has a default implementation of value which returns zero, as time-dependent simulation based optimization problems may have an objective value which depends only on the final state of the system.
void advanceAdjoint(Vector< Real > &l, Vector< Real > &rhs, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
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 solveState(const Vector< Real > &x)
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Real updateSketch(const Vector< Real > &x, const Real tol)
typename std::vector< Real >::size_type size_type
PartitionedVector< Real > & partition(Vector< Real > &x) const
Ptr< Vector< Real > > makeDynamicVector(const Vector< Real > &x) const
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)
const Ptr< DynamicObjective< Real > > obj_
const Ptr< DynamicConstraint< Real > > con_
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 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)
std::vector< Ptr< Vector< Real > > > whist_
void solveAdjoint(const Vector< Real > &x)
void setTerminalCondition(Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)
std::vector< Ptr< Vector< Real > > > uhist_
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)
std::vector< Ptr< Vector< Real > > > lhist_
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)
std::vector< Ptr< Vector< Real > > > phist_
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
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)
ROL::Ptr< const Vector< Real > > get(size_type i) const
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)
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
const PartitionedVector< Real > & partition(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 updateGradient(Vector< Real > &g, const Vector< Real > &l, const Vector< Real > &uold, const Vector< Real > &unew, const Vector< Real > &z, const TimeStamp< Real > &ts)