ROL
ROL_Reduced_Objective_SimOpt.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
44 
45 #ifndef ROL_REDUCED_OBJECTIVE_SIMOPT_H
46 #define ROL_REDUCED_OBJECTIVE_SIMOPT_H
47 
48 #include "ROL_Objective_SimOpt.hpp"
50 
51 namespace ROL {
52 
53 template <class Real>
54 class Reduced_Objective_SimOpt : public Objective<Real> {
55 private:
56  Teuchos::RCP<Objective_SimOpt<Real> > obj_;
57  Teuchos::RCP<EqualityConstraint_SimOpt<Real> > con_;
58 
59  // Primal vectors
60  Teuchos::RCP<Vector<Real> > state_;
61  Teuchos::RCP<Vector<Real> > state_sens_;
62  Teuchos::RCP<Vector<Real> > adjoint_;
63  Teuchos::RCP<Vector<Real> > adjoint_sens_;
64 
65  // Dual vectors
66  Teuchos::RCP<Vector<Real> > dualstate_;
67  Teuchos::RCP<Vector<Real> > dualstate1_;
68  Teuchos::RCP<Vector<Real> > dualadjoint_;
69  Teuchos::RCP<Vector<Real> > dualcontrol_;
70 
71  // Vector storage
72  std::map<std::vector<Real>,Teuchos::RCP<Vector<Real> > > state_storage_;
73  std::map<std::vector<Real>,Teuchos::RCP<Vector<Real> > > adjoint_storage_;
74 
75  bool storage_;
78 
79  void solve_state_equation(const Vector<Real> &x, Real &tol, bool flag = true, int iter = -1) {
80  // Solve state equation if not done already
81  if ( state_storage_.count(this->getParameter()) && storage_ ) {
82  state_->set(*state_storage_[this->getParameter()]);
83  }
84  else {
85  // Update equality constraint
86  con_->update_2(x,flag,iter);
87  // Solve state
88  con_->solve(*dualadjoint_,*state_,x,tol);
89  // Update equality constraint
90  con_->update_1(*state_,flag,iter);
91  // Update full objective function
92  obj_->update(*state_,x,flag,iter);
93  // Store state
94  if ( storage_ ) {
95  state_storage_.insert(
96  std::pair<std::vector<Real>,Teuchos::RCP<Vector<Real> > >(
97  this->getParameter(),state_->clone()));
98  state_storage_[this->getParameter()]->set(*state_);
99  }
100  }
101  }
102 
107  void solve_adjoint_equation(const Vector<Real> &x, Real &tol) {
108  // Solve adjoint equation if not done already
109  if ( adjoint_storage_.count(this->getParameter()) && storage_ ) {
110  adjoint_->set(*adjoint_storage_[this->getParameter()]);
111  }
112  else {
113  // Evaluate the full gradient wrt u
114  obj_->gradient_1(*dualstate_,*state_,x,tol);
115  // Solve adjoint equation
116  con_->applyInverseAdjointJacobian_1(*adjoint_,*dualstate_,*state_,x,tol);
117  adjoint_->scale(-1.0);
118  // Store adjoint
119  if ( storage_ ) {
120  adjoint_storage_.insert(
121  std::pair<std::vector<Real>,Teuchos::RCP<Vector<Real> > >(
122  this->getParameter(),adjoint_->clone()));
123  adjoint_storage_[this->getParameter()]->set(*adjoint_);
124  }
125  }
126  }
127 
132  void solve_state_sensitivity(const Vector<Real> &v, const Vector<Real> &x, Real &tol) {
133  // Solve state sensitivity equation
134  con_->applyJacobian_2(*dualadjoint_,v,*state_,x,tol);
135  dualadjoint_->scale(-1.0);
136  con_->applyInverseJacobian_1(*state_sens_,*dualadjoint_,*state_,x,tol);
137  }
138 
146  void solve_adjoint_sensitivity(const Vector<Real> &v, const Vector<Real> &x, Real &tol) {
147  // Evaluate full hessVec in the direction (s,v)
148  obj_->hessVec_11(*dualstate_,*state_sens_,*state_,x,tol);
149  obj_->hessVec_12(*dualstate1_,v,*state_,x,tol);
150  dualstate_->plus(*dualstate1_);
151  // Apply adjoint Hessian of constraint
152  con_->applyAdjointHessian_11(*dualstate1_,*adjoint_,*state_sens_,*state_,x,tol);
153  dualstate_->plus(*dualstate1_);
154  con_->applyAdjointHessian_21(*dualstate1_,*adjoint_,v,*state_,x,tol);
155  dualstate_->plus(*dualstate1_);
156  // Solve adjoint sensitivity equation
157  dualstate_->scale(-1.0);
158  con_->applyInverseAdjointJacobian_1(*adjoint_sens_,*dualstate_,*state_,x,tol);
159  }
160 
161 public:
171  const Teuchos::RCP<EqualityConstraint_SimOpt<Real> > &con,
172  const Teuchos::RCP<Vector<Real> > &state,
173  const Teuchos::RCP<Vector<Real> > &adjoint,
174  const bool storage = true,
175  const bool useFDhessVec = false)
176  : obj_(obj), con_(con),
177  state_(state), state_sens_(state->clone()),
178  adjoint_(adjoint), adjoint_sens_(adjoint->clone()),
179  dualstate_(state->dual().clone()), dualstate1_(state->dual().clone()),
180  dualadjoint_(adjoint->dual().clone()), dualcontrol_(Teuchos::null),
181  storage_(storage), useFDhessVec_(useFDhessVec), is_initialized_(false) {
182  state_storage_.clear();
183  adjoint_storage_.clear();
184  }
185 
199  const Teuchos::RCP<EqualityConstraint_SimOpt<Real> > &con,
200  const Teuchos::RCP<Vector<Real> > &state,
201  const Teuchos::RCP<Vector<Real> > &adjoint,
202  const Teuchos::RCP<Vector<Real> > &dualstate,
203  const Teuchos::RCP<Vector<Real> > &dualadjoint,
204  const bool storage = true,
205  const bool useFDhessVec = false)
206  : obj_(obj), con_(con),
207  state_(state), state_sens_(state->clone()),
208  adjoint_(adjoint), adjoint_sens_(adjoint->clone()),
209  dualstate_(dualstate), dualstate1_(dualstate->clone()),
210  dualadjoint_(dualadjoint->clone()), dualcontrol_(Teuchos::null),
211  storage_(storage), useFDhessVec_(useFDhessVec), is_initialized_(false) {
212  state_storage_.clear();
213  adjoint_storage_.clear();
214  }
215 
218  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
219  // Reset storage flags
220  state_storage_.clear();
221  if ( flag ) {
222  adjoint_storage_.clear();
223  }
224  }
225 
230  Real value( const Vector<Real> &x, Real &tol ) {
231  // Solve state equation
232  solve_state_equation(x,tol);
233  // Get objective function value
234  return obj_->value(*state_,x,tol);
235  }
236 
242  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
243  if (!is_initialized_) {
244  dualcontrol_ = g.clone();
245  is_initialized_ = true;
246  }
247  // Solve state equation
248  solve_state_equation(x,tol);
249  // Solve adjoint equation
250  solve_adjoint_equation(x,tol);
251  // Evaluate the full gradient wrt z
252  obj_->gradient_2(*dualcontrol_,*state_,x,tol);
253  // Build gradient
254  con_->applyAdjointJacobian_2(g,*adjoint_,*state_,x,tol);
255  g.plus(*dualcontrol_);
256  }
257 
261  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
262  if ( useFDhessVec_ ) {
263  Objective<Real>::hessVec(hv,v,x,tol);
264  }
265  else {
266  if (!is_initialized_) {
267  dualcontrol_ = hv.clone();
268  is_initialized_ = true;
269  }
270  // Solve state equation
271  solve_state_equation(x,tol);
272  // Solve adjoint equation
273  solve_adjoint_equation(x,tol);
274  // Solve state sensitivity equation
275  solve_state_sensitivity(v,x,tol);
276  // Solve adjoint sensitivity equation
277  solve_adjoint_sensitivity(v,x,tol);
278  // Build hessVec
279  con_->applyAdjointJacobian_2(hv,*adjoint_sens_,*state_,x,tol);
280  obj_->hessVec_21(*dualcontrol_,*state_sens_,*state_,x,tol);
281  hv.plus(*dualcontrol_);
282  obj_->hessVec_22(*dualcontrol_,v,*state_,x,tol);
283  hv.plus(*dualcontrol_);
284  con_->applyAdjointHessian_12(*dualcontrol_,*adjoint_,*state_sens_,*state_,x,tol);
285  hv.plus(*dualcontrol_);
286  con_->applyAdjointHessian_22(*dualcontrol_,*adjoint_,v,*state_,x,tol);
287  hv.plus(*dualcontrol_);
288  }
289  }
290 
293  virtual void precond( Vector<Real> &Pv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
294  Pv.set(v.dual());
295  }
296 
297 // For parametrized (stochastic) objective functions and constraints
298 public:
299  void setParameter(const std::vector<Real> &param) {
301  con_->setParameter(param);
302  obj_->setParameter(param);
303  }
304 }; // class Reduced_Objective_SimOpt
305 
306 } // namespace ROL
307 
308 #endif
Teuchos::RCP< Objective_SimOpt< Real > > obj_
Provides the interface to evaluate objective functions.
Provides the interface to evaluate simulation-based objective functions.
virtual void plus(const Vector &x)=0
Compute , where .
virtual void precond(Vector< Real > &Pv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply a reduced Hessian preconditioner.
void solve_state_equation(const Vector< Real > &x, Real &tol, bool flag=true, int iter=-1)
Reduced_Objective_SimOpt(const Teuchos::RCP< Objective_SimOpt< Real > > &obj, const Teuchos::RCP< EqualityConstraint_SimOpt< Real > > &con, const Teuchos::RCP< Vector< Real > > &state, const Teuchos::RCP< Vector< Real > > &adjoint, const bool storage=true, const bool useFDhessVec=false)
Constructor.
void setParameter(const std::vector< Real > &param)
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Teuchos::RCP< Vector< Real > > state_
Teuchos::RCP< Vector< Real > > adjoint_sens_
Teuchos::RCP< Vector< Real > > state_sens_
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Given , evaluate the gradient of the objective function where solves .
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Defines the equality constraint operator interface for simulation-based optimization.
Teuchos::RCP< Vector< Real > > adjoint_
Teuchos::RCP< Vector< Real > > dualstate_
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:213
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Given , evaluate the Hessian of the objective function in the direction .
std::map< std::vector< Real >, Teuchos::RCP< Vector< Real > > > adjoint_storage_
std::map< std::vector< Real >, Teuchos::RCP< Vector< Real > > > state_storage_
Real value(const Vector< Real > &x, Real &tol)
Given , evaluate the objective function where solves .
const std::vector< Real > getParameter(void) const
virtual void setParameter(const std::vector< Real > &param)
Teuchos::RCP< Vector< Real > > dualcontrol_
void solve_adjoint_sensitivity(const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Given , the adjoint variable , and a direction , solve the adjoint sensitvity equation for ...
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
void solve_state_sensitivity(const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Given which solves the state equation and a direction , solve the state senstivity equation for ...
Teuchos::RCP< EqualityConstraint_SimOpt< Real > > con_
Teuchos::RCP< Vector< Real > > dualstate1_
Teuchos::RCP< Vector< Real > > dualadjoint_
Reduced_Objective_SimOpt(const Teuchos::RCP< Objective_SimOpt< Real > > &obj, const Teuchos::RCP< EqualityConstraint_SimOpt< Real > > &con, const Teuchos::RCP< Vector< Real > > &state, const Teuchos::RCP< Vector< Real > > &adjoint, const Teuchos::RCP< Vector< Real > > &dualstate, const Teuchos::RCP< Vector< Real > > &dualadjoint, const bool storage=true, const bool useFDhessVec=false)
Secondary, general constructor for use with dual optimization vector spaces where the user does not d...
void solve_adjoint_equation(const Vector< Real > &x, Real &tol)
Given which solves the state equation, solve the adjoint equation for .
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update the SimOpt objective function and equality constraint.