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  gp.get(0)->zero(); // zero for the nonexistant zeroth control interval
252  const PartitionedVector<Real> &xp = partition(x);
253  const Real one(1);
254  size_type uindex = (useSketch_ ? 1 : Nt_-1);
255  size_type lindex = (useSketch_ ? 0 : Nt_-1);
256  // Must first compute the value
257  solveState(x);
258  if (useSketch_ && useInexact_) {
259  tol = updateSketch(x,tol);
260  }
261  // Recover state from sketch
262  if (useSketch_) {
263  uhist_[1]->set(*uhist_[0]);
264  stateSketch_->reconstruct(*uhist_[0],static_cast<int>(Nt_)-3);
265  if (isAdjointComputed_) {
266  adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(Nt_)-2);
267  }
268  }
269  // Update dynamic constraint and objective
270  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
271  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
272  // Solve for terminal condition
273  if (!isAdjointComputed_) {
274  setTerminalCondition(*lhist_[lindex],
275  *uhist_[uindex-1], *uhist_[uindex],
276  *xp.get(Nt_-1), timeStamp_[Nt_-1]);
277  if (useSketch_) {
278  adjointSketch_->advance(one,*lhist_[0],static_cast<int>(Nt_)-2,one);
279  }
280  }
281  // Update gradient on terminal interval
282  updateGradient(*gp.get(Nt_-1), *lhist_[lindex],
283  *uhist_[uindex-1], *uhist_[uindex],
284  *xp.get(Nt_-1), timeStamp_[Nt_-1]);
285  // Run reverse time stepper
286  for (size_type k = Nt_-2; k > 0; --k) {
287  if (!isAdjointComputed_) {
288  // Compute k+1 component of rhs
289  computeAdjointRHS(*rhs_, *lhist_[lindex],
290  *uhist_[uindex-1], *uhist_[uindex],
291  *xp.get(k+1), timeStamp_[k+1]);
292  }
293  uindex = (useSketch_ ? 1 : k);
294  lindex = (useSketch_ ? 0 : k);
295  // Recover state from sketch
296  if (useSketch_) {
297  uhist_[1]->set(*uhist_[0]);
298  if (k==1) {
299  uhist_[0]->set(*u0_);
300  }
301  else {
302  stateSketch_->reconstruct(*uhist_[0],static_cast<int>(k)-2);
303  }
304  if (isAdjointComputed_) {
305  adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(k)-1);
306  }
307  }
308  // Update dynamic constraint and objective
309  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
310  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
311  // Solve for adjoint on interval k
312  if (!isAdjointComputed_) {
313  advanceAdjoint(*lhist_[lindex], *rhs_,
314  *uhist_[uindex-1], *uhist_[uindex],
315  *xp.get(k), timeStamp_[k]);
316  if (useSketch_) {
317  adjointSketch_->advance(one,*lhist_[0],static_cast<int>(k)-1,one);
318  }
319  }
320  // Update gradient on interval k
321  updateGradient(*gp.get(k), *lhist_[lindex],
322  *uhist_[uindex-1], *uhist_[uindex],
323  *xp.get(k), timeStamp_[k]);
324  }
325  isAdjointComputed_ = true;
326  }
327 
328  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
329  if (useHessian_) {
330  // Must first solve the state and adjoint equations
331  solveState(x);
332  solveAdjoint(x);
333  // Now compute Hessian
334  const Real one(1);
335  const PartitionedVector<Real> &xp = partition(x);
336  const PartitionedVector<Real> &vp = partition(v);
338  hvp.get(0)->zero(); // zero for the nonexistant zeroth control interval
339  // Compute state sensitivity
340  whist_[0]->zero();
341  if (useSketch_) {
342  stateSensSketch_->update();
343  //stateSensSketch_->advance(one,*whist_[0],0,one);
344  uhist_[0]->set(*u0_);
345  }
346  size_type uindex, lindex;
347  for (size_type k = 1; k < Nt_; ++k) {
348  uindex = (useSketch_ ? 1 : k);
349  lindex = (useSketch_ ? 0 : k);
350  // Reconstruct sketched state
351  if (useSketch_) {
352  stateSketch_->reconstruct(*uhist_[1],static_cast<int>(k)-1);
353  adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(k)-1);
354  }
355  // Update dynamic constraint and objective
356  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
357  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
358  // Advance state sensitivity on current time interval
359  advanceStateSens(*whist_[uindex],
360  *vp.get(k), *whist_[uindex-1],
361  *uhist_[uindex-1], *uhist_[uindex],
362  *xp.get(k), timeStamp_[k]);
363  // Set control only Hessian
364  computeControlHessLag(*hvp.get(k),
365  *vp.get(k), *lhist_[lindex],
366  *uhist_[uindex-1], *uhist_[uindex],
367  *xp.get(k), timeStamp_[k]);
368  if (!useSymHess_) {
369  // Add mixed derivative Hessian
370  addMixedHessLag(*hvp.get(k), *lhist_[lindex],
371  *whist_[uindex-1], *whist_[uindex],
372  *uhist_[uindex-1], *uhist_[uindex],
373  *xp.get(k), timeStamp_[k]);
374  }
375  // Sketch state sensitivity
376  if (useSketch_) {
377  stateSensSketch_->advance(one,*whist_[1],static_cast<int>(k)-1,one);
378  whist_[0]->set(*whist_[1]);
379  uhist_[0]->set(*uhist_[1]);
380  }
381  }
382 
383  // Compute adjoint sensitivity
384  uindex = (useSketch_ ? 1 : Nt_-1);
385  lindex = (useSketch_ ? 0 : Nt_-1);
386  if (useSketch_) {
387  // Recover terminal state
388  uhist_[1]->set(*uhist_[0]);
389  stateSketch_->reconstruct(*uhist_[0],static_cast<int>(Nt_)-3);
390  // Recover terminal adjoint
391  adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(Nt_)-2);
392  // Recover terminal state sensitivity
393  whist_[1]->set(*whist_[0]);
394  stateSensSketch_->reconstruct(*whist_[0],static_cast<int>(Nt_)-3);
395  }
396  // Update dynamic constraint and objective
397  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
398  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
399  // Solve for terminal condition
401  *vp.get(Nt_-1), *lhist_[lindex],
402  *whist_[uindex-1], *whist_[uindex],
403  *uhist_[uindex-1], *uhist_[uindex],
404  *xp.get(Nt_-1), timeStamp_[Nt_-1]);
405  if (useSymHess_) {
406  // Add mixed derivative Hessian
407  addMixedHessLag(*hvp.get(Nt_-1), *lhist_[lindex],
408  *whist_[uindex-1], *whist_[uindex],
409  *uhist_[uindex-1], *uhist_[uindex],
410  *xp.get(Nt_-1), timeStamp_[Nt_-1]);
411  }
412  // Add adjoint sensitivity to Hessian
413  addAdjointSens(*hvp.get(Nt_-1), *phist_[lindex],
414  *uhist_[uindex-1], *uhist_[uindex],
415  *xp.get(Nt_-1), timeStamp_[Nt_-1]);
416  // Run reverse time stepper
417  for (size_type k = Nt_-2; k > 0; --k) {
418  // Compute new components of rhs
420  *whist_[uindex-1], *whist_[uindex],
421  *uhist_[uindex-1], *uhist_[uindex],
422  *xp.get(k+1), timeStamp_[k+1],
423  false);
425  *vp.get(k+1), *lhist_[lindex],
426  *uhist_[uindex-1], *uhist_[uindex],
427  *xp.get(k+1), timeStamp_[k+1],
428  true);
430  *uhist_[uindex-1], *uhist_[uindex],
431  *xp.get(k+1), timeStamp_[k+1],
432  true);
433  // Recover state, adjoint and state sensitivity from sketch
434  if (useSketch_) {
435  uhist_[1]->set(*uhist_[0]);
436  whist_[1]->set(*whist_[0]);
437  if (k==1) {
438  uhist_[0]->set(*u0_);
439  whist_[0]->zero();
440  }
441  else {
442  stateSketch_->reconstruct(*uhist_[0],static_cast<int>(k)-2);
443  stateSensSketch_->reconstruct(*whist_[0],static_cast<int>(k)-2);
444  }
445  adjointSketch_->reconstruct(*lhist_[0],static_cast<int>(k)-1);
446  }
447  uindex = (useSketch_ ? 1 : k);
448  lindex = (useSketch_ ? 0 : k);
449  // Update dynamic constraint and objective
450  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
451  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
452  // Compute old components of rhs
454  *whist_[uindex-1], *whist_[uindex],
455  *uhist_[uindex-1], *uhist_[uindex],
456  *xp.get(k), timeStamp_[k],
457  true);
459  *vp.get(k), *lhist_[lindex],
460  *uhist_[uindex-1], *uhist_[uindex],
461  *xp.get(k), timeStamp_[k],
462  true);
463  // Solve for adjoint on interval k
464  advanceAdjointSens(*phist_[lindex], *rhs_,
465  *uhist_[uindex-1], *uhist_[uindex],
466  *xp.get(k), timeStamp_[k]);
467  if (useSymHess_) {
468  // Add mixed derivative Hessian
469  addMixedHessLag(*hvp.get(k), *lhist_[lindex],
470  *whist_[uindex-1], *whist_[uindex],
471  *uhist_[uindex-1], *uhist_[uindex],
472  *xp.get(k), timeStamp_[k]);
473  }
474  // Add adjoint sensitivity to Hessian
475  addAdjointSens(*hvp.get(k), *phist_[lindex],
476  *uhist_[uindex-1], *uhist_[uindex],
477  *xp.get(k), timeStamp_[k]);
478  }
479  }
480  else {
481  Objective<Real>::hessVec(hv,v,x,tol);
482  }
483  }
484 
485 private:
486  /***************************************************************************/
487  /************ Method to solve state equation *******************************/
488  /***************************************************************************/
489  void solveState(const Vector<Real> &x) {
490  if (!isStateComputed_) {
491  const Real one(1);
492  const PartitionedVector<Real> &xp = partition(x);
493  // Set initial condition
494  uhist_[0]->set(*u0_);
495  if (useSketch_) {
496  //stateSketch_->advance(one,*uhist_[0],0,one);
497  }
498  // Run time stepper
499  size_type index;
500  for (size_type k = 1; k < Nt_; ++k) {
501  index = (useSketch_ ? 1 : k);
502  // Update dynamic constraint and objective
503  con_->update(*uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
504  obj_->update(*uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
505  // Solve state on current time interval
506  con_->solve(*cprimal_, *uhist_[index-1], *uhist_[index], *xp.get(k), timeStamp_[k]);
507  // Sketch state
508  if (useSketch_) {
509  stateSketch_->advance(one,*uhist_[index],static_cast<int>(k)-1,one);
510  uhist_[index-1]->set(*uhist_[index]);
511  }
512  }
513  isStateComputed_ = true;
514  }
515  }
516 
517  Real updateSketch(const Vector<Real> &x, const Real tol) {
518  const PartitionedVector<Real> &xp = partition(x);
519  Real err(0), cnorm(0);
520  bool flag = true;
521  while (flag) {
522  err = static_cast<Real>(0);
523  uhist_[0]->set(*u0_);
524  for (size_type k = 1; k < Nt_; ++k) {
525  stateSketch_->reconstruct(*uhist_[1],static_cast<int>(k));
526  con_->value(*cprimal_, *uhist_[0], *uhist_[1], *xp.get(k), timeStamp_[k]);
527  cnorm = cprimal_->norm();
528  err = (cnorm > err ? cnorm : err); // Use Linf norm. Could change to use, e.g., L2.
529  if (err > tol) {
530  break;
531  }
532  }
533  if (err > tol) {
534  rankState_ *= updateFactor_; // Perhaps there is a better update strategy
536  stateSketch_->setRank(rankState_);
537  solveState(x);
538  }
539  else {
540  flag = false;
541  break;
542  }
543  }
544  return err;
545  }
546 
547  /***************************************************************************/
548  /************ Methods to solve adjoint equation ****************************/
549  /***************************************************************************/
551  const Vector<Real> &uold, const Vector<Real> &unew,
552  const Vector<Real> &z, const TimeStamp<Real> &ts) {
553  obj_->gradient_un(*udual_, uold, unew, z, ts);
554  con_->applyInverseAdjointJacobian_un(l, *udual_, uold, unew, z, ts);
555  }
556 
558  const Vector<Real> &uold, const Vector<Real> &unew,
559  const Vector<Real> &z, const TimeStamp<Real> &ts) {
560  const Real one(1);
561  obj_->gradient_uo(rhs, uold, unew, z, ts);
562  con_->applyAdjointJacobian_uo(*udual_, l, uold, unew, z, ts);
563  rhs.axpy(-one,*udual_);
564  }
565 
567  const Vector<Real> &uold, const Vector<Real> &unew,
568  const Vector<Real> &z, const TimeStamp<Real> &ts) {
569  obj_->gradient_un(*udual_, uold, unew, z, ts);
570  rhs.plus(*udual_);
571  con_->applyInverseAdjointJacobian_un(l, rhs, uold, unew, z, ts);
572  }
573 
574  void solveAdjoint(const Vector<Real> &x) {
575  if (!isAdjointComputed_) {
576  const PartitionedVector<Real> &xp = partition(x);
577  const Real one(1);
578  size_type uindex = (useSketch_ ? 1 : Nt_-1);
579  size_type lindex = (useSketch_ ? 0 : Nt_-1);
580  // Must first compute solve the state equation
581  solveState(x);
582  // Recover state from sketch
583  if (useSketch_) {
584  uhist_[1]->set(*uhist_[0]);
585  stateSketch_->reconstruct(*uhist_[0],static_cast<int>(Nt_)-3);
586  }
587  // Update dynamic constraint and objective
588  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
589  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(Nt_-1), timeStamp_[Nt_-1]);
590  // Solve for terminal condition
591  setTerminalCondition(*lhist_[lindex],
592  *uhist_[uindex-1], *uhist_[uindex],
593  *xp.get(Nt_-1), timeStamp_[Nt_-1]);
594  if (useSketch_) {
595  adjointSketch_->advance(one,*lhist_[lindex],static_cast<int>(Nt_)-2,one);
596  }
597  // Run reverse time stepper
598  for (size_type k = Nt_-2; k > 0; --k) {
599  // Compute k+1 component of rhs
600  computeAdjointRHS(*rhs_, *lhist_[lindex],
601  *uhist_[uindex-1], *uhist_[uindex],
602  *xp.get(k+1), timeStamp_[k+1]);
603  // Recover state from sketch
604  if (useSketch_) {
605  uhist_[1]->set(*uhist_[0]);
606  if (k==1) {
607  uhist_[0]->set(*u0_);
608  }
609  else {
610  stateSketch_->reconstruct(*uhist_[0],static_cast<int>(k)-2);
611  }
612  }
613  uindex = (useSketch_ ? 1 : k);
614  lindex = (useSketch_ ? 0 : k);
615  // Update dynamic constraint and objective
616  con_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
617  obj_->update(*uhist_[uindex-1], *uhist_[uindex], *xp.get(k), timeStamp_[k]);
618  // Solve for adjoint on interval k
619  advanceAdjoint(*lhist_[lindex], *rhs_,
620  *uhist_[uindex-1], *uhist_[uindex],
621  *xp.get(k), timeStamp_[k]);
622  if (useSketch_) {
623  adjointSketch_->advance(one,*lhist_[lindex],static_cast<int>(k)-1,one);
624  }
625  }
626  isAdjointComputed_ = true;
627  }
628  }
629 
630  /***************************************************************************/
631  /************ Method for gradient computation ******************************/
632  /***************************************************************************/
634  const Vector<Real> &uold, const Vector<Real> &unew,
635  const Vector<Real> &z, const TimeStamp<Real> &ts) {
636  const Real one(1);
637  obj_->gradient_z(g, uold, unew, z, ts);
638  con_->applyAdjointJacobian_z(*zdual_, l, uold, unew, z, ts);
639  g.axpy(-one,*zdual_);
640  }
641 
642  /***************************************************************************/
643  /************ Method to solve state sensitivity equation *******************/
644  /***************************************************************************/
646  const Vector<Real> &wold, const Vector<Real> &uold,
647  const Vector<Real> &unew, const Vector<Real> &z,
648  const TimeStamp<Real> &ts) {
649  const Real one(1);
650  con_->applyJacobian_z(*crhs_, v, uold, unew, z, ts);
651  con_->applyJacobian_uo(*cprimal_, wold, uold, unew, z, ts);
652  crhs_->axpy(-one, *cprimal_);
653  con_->applyInverseJacobian_un(wnew, *crhs_, uold, unew, z, ts);
654  }
655 
656  /***************************************************************************/
657  /************ Methods to solve adjoint sensitivity equation ****************/
658  /***************************************************************************/
660  const Vector<Real> &v, const Vector<Real> &l,
661  const Vector<Real> &wold, const Vector<Real> &wnew,
662  const Vector<Real> &uold, const Vector<Real> &unew,
663  const Vector<Real> &z, const TimeStamp<Real> &ts) {
664  const Real one(1);
665  // Mixed derivative rhs term
666  con_->applyAdjointHessian_z_un(*rhs_, l, v, uold, unew, z, ts);
667  obj_->hessVec_un_z(*udual_, v, uold, unew, z, ts);
668  rhs_->axpy(-one, *udual_);
669  // State derivative rhs term
670  con_->applyAdjointHessian_un_un(*udual_, l, wnew, uold, unew, z, ts);
671  rhs_->axpy(-one, *udual_);
672  obj_->hessVec_un_un(*udual_, wnew, uold, unew, z, ts);
673  rhs_->plus(*udual_);
674  con_->applyAdjointHessian_uo_un(*udual_, l, wold, uold, unew, z, ts);
675  rhs_->axpy(-one, *udual_);
676  obj_->hessVec_un_uo(*udual_, wold, uold, unew, z, ts);
677  rhs_->plus(*udual_);
678  // Invert adjoint Jacobian
679  con_->applyInverseAdjointJacobian_un(p, *rhs_, uold, unew, z, ts);
680  }
681 
683  const Vector<Real> &wold, const Vector<Real> &wnew,
684  const Vector<Real> &uold, const Vector<Real> &unew,
685  const Vector<Real> &z, const TimeStamp<Real> &ts,
686  const bool sumInto = false) {
687  const Real one(1);
688  // Compute new/old Hessian of Lagrangian
689  if (!sumInto) {
690  obj_->hessVec_un_uo(Hv, wold, uold, unew, z, ts);
691  }
692  else {
693  obj_->hessVec_un_uo(*udual_, wold, uold, unew, z, ts);
694  Hv.plus(*udual_);
695  }
696  con_->applyAdjointHessian_uo_un(*udual_, l, wold, uold, unew, z, ts);
697  Hv.axpy(-one,*udual_);
698  // Compute new/new Hessian of Lagrangian
699  obj_->hessVec_un_un(*udual_, wnew, uold, unew, z, ts);
700  Hv.plus(*udual_);
701  con_->applyAdjointHessian_un_un(*udual_, l, wnew, uold, unew, z, ts);
702  Hv.axpy(-one,*udual_);
703  }
704 
706  const Vector<Real> &v, const Vector<Real> &l,
707  const Vector<Real> &uold, const Vector<Real> &unew,
708  const Vector<Real> &z, const TimeStamp<Real> &ts,
709  const bool sumInto = false) {
710  const Real one(1);
711  // Compute new/old Hessian of Lagrangian
712  if (!sumInto) {
713  con_->applyAdjointHessian_z_un(Hv, l, v, uold, unew, z, ts);
714  }
715  else {
716  con_->applyAdjointHessian_z_un(*udual_, l, v, uold, unew, z, ts);
717  Hv.plus(*udual_);
718  }
719  obj_->hessVec_un_z(*udual_, v, uold, unew, z, ts);
720  Hv.axpy(-one, *udual_);
721  }
722 
724  const Vector<Real> &uold, const Vector<Real> &unew,
725  const Vector<Real> &z, const TimeStamp<Real> &ts,
726  const bool sumInto = false) {
727  const Real one(1);
728  if (!sumInto) {
729  con_->applyAdjointJacobian_uo(Hv, p, uold, unew, z, ts);
730  Hv.scale(-one);
731  }
732  else {
733  con_->applyAdjointJacobian_uo(*udual_, p, uold, unew, z, ts);
734  Hv.axpy(-one, *udual_);
735  }
736  }
737 
739  const Vector<Real> &wold, const Vector<Real> &wnew,
740  const Vector<Real> &uold, const Vector<Real> &unew,
741  const Vector<Real> &z, const TimeStamp<Real> &ts,
742  const bool sumInto = false) {
743  const Real one(1);
744  // Compute old/new Hessian of Lagrangian
745  if (!sumInto) {
746  obj_->hessVec_uo_un(Hv, wnew, uold, unew, z, ts);
747  }
748  else {
749  obj_->hessVec_uo_un(*udual_, wnew, uold, unew, z, ts);
750  Hv.plus(*udual_);
751  }
752  con_->applyAdjointHessian_un_uo(*udual_, l, wnew, uold, unew, z, ts);
753  Hv.axpy(-one,*udual_);
754  // Compute old/old Hessian of Lagrangian
755  obj_->hessVec_uo_uo(*udual_, wold, uold, unew, z, ts);
756  Hv.plus(*udual_);
757  con_->applyAdjointHessian_uo_uo(*udual_, l, wold, uold, unew, z, ts);
758  Hv.axpy(-one,*udual_);
759  }
760 
762  const Vector<Real> &v, const Vector<Real> &l,
763  const Vector<Real> &uold, const Vector<Real> &unew,
764  const Vector<Real> &z, const TimeStamp<Real> &ts,
765  const bool sumInto = false) {
766  const Real one(1);
767  // Compute new/old Hessian of Lagrangian
768  if (!sumInto) {
769  con_->applyAdjointHessian_z_uo(Hv, l, v, uold, unew, z, ts);
770  }
771  else {
772  con_->applyAdjointHessian_z_uo(*udual_, l, v, uold, unew, z, ts);
773  Hv.plus(*udual_);
774  }
775  obj_->hessVec_uo_z(*udual_, v, uold, unew, z, ts);
776  Hv.axpy(-one,*udual_);
777  }
778 
780  const Vector<Real> &uold, const Vector<Real> &unew,
781  const Vector<Real> &z, const TimeStamp<Real> &ts) {
782  // Solve adjoint sensitivity on current time interval
783  con_->applyInverseAdjointJacobian_un(p, rhs, uold, unew, z, ts);
784  }
785 
786  /***************************************************************************/
787  /************ Method for Hessian-times-a-vector computation ****************/
788  /***************************************************************************/
790  const Vector<Real> &v, const Vector<Real> &l,
791  const Vector<Real> &uold, const Vector<Real> &unew,
792  const Vector<Real> &z, const TimeStamp<Real> &ts) {
793  const Real one(1);
794  // Compute Hessian of Lagrangian
795  obj_->hessVec_z_z(Hv, v, uold, unew, z, ts);
796  con_->applyAdjointHessian_z_z(*zdual_, l, v, uold, unew, z, ts);
797  Hv.axpy(-one, *zdual_);
798  }
799 
801  const Vector<Real> &wold, const Vector<Real> &wnew,
802  const Vector<Real> &uold, const Vector<Real> &unew,
803  const Vector<Real> &z, const TimeStamp<Real> &ts) {
804  const Real one(1);
805  // Compute Hessian of Lagrangian on previous time interval
806  obj_->hessVec_z_uo(*zdual_, wold, uold, unew, z, ts);
807  Hv.axpy(-one, *zdual_);
808  con_->applyAdjointHessian_uo_z(*zdual_, l, wold, uold, unew, z, ts);
809  Hv.plus(*zdual_);
810  // Compute Hessian of Lagrangian on current time interval
811  obj_->hessVec_z_un(*zdual_, wnew, uold, unew, z, ts);
812  Hv.axpy(-one, *zdual_);
813  con_->applyAdjointHessian_un_z(*zdual_, l, wnew, uold, unew, z, ts);
814  Hv.plus(*zdual_);
815  }
816 
818  const Vector<Real> &uold, const Vector<Real> &unew,
819  const Vector<Real> &z, const TimeStamp<Real> &ts) {
820  con_->applyAdjointJacobian_z(*zdual_, p, uold, unew, z, ts);
821  Hv.plus(*zdual_);
822  }
823 };
824 
825 } // namespace ROL
826 
827 #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)
const PartitionedVector< Real > & partition(const Vector< Real > &x) const
ROL::Ptr< const Vector< Real > > get(size_type i) const
Defines the reduced time-dependent objective function interface for simulation-based optimization...
virtual void plus(const Vector &x)=0
Compute , where .
Ptr< Vector< Real > > makeDynamicVector(const Vector< Real > &x) const
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.
PartitionedVector< Real > & partition(Vector< Real > &x) const
Real updateSketch(const Vector< Real > &x, const Real tol)
typename std::vector< Real >::size_type size_type
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)
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.
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)