ROL
ROL_Reduced_Objective_SimOpt_Def.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_REDUCED_OBJECTIVE_SIMOPT_DEF_H
45 #define ROL_REDUCED_OBJECTIVE_SIMOPT_DEF_H
46 
47 namespace ROL {
48 
49 template<typename Real>
51  const Ptr<Objective_SimOpt<Real>> &obj,
52  const Ptr<Constraint_SimOpt<Real>> &con,
53  const Ptr<Vector<Real>> &state,
54  const Ptr<Vector<Real>> &control,
55  const Ptr<Vector<Real>> &adjoint,
56  const bool storage,
57  const bool useFDhessVec)
58  : obj_(obj), con_(con),
59  storage_(storage), useFDhessVec_(useFDhessVec),
60  nupda_(0), nvalu_(0), ngrad_(0), nhess_(0), nprec_(0),
61  nstat_(0), nadjo_(0), nssen_(0), nasen_(0),
62  updateFlag_(true), updateIter_(0), updateType_(UpdateType::Initial),
63  newUpdate_(false) {
64  stateStore_ = makePtr<VectorController<Real>>();
65  adjointStore_ = makePtr<VectorController<Real>>();
66  state_ = state->clone(); state_->set(*state);
67  adjoint_ = adjoint->clone();
68  state_sens_ = state->clone();
69  adjoint_sens_ = adjoint->clone();
70  dualstate_ = state->dual().clone();
71  dualstate1_ = state->dual().clone();
72  dualadjoint_ = adjoint->dual().clone();
73  dualcontrol_ = control->dual().clone();
74 }
75 
76 template<typename Real>
78  const Ptr<Objective_SimOpt<Real>> &obj,
79  const Ptr<Constraint_SimOpt<Real>> &con,
80  const Ptr<Vector<Real>> &state,
81  const Ptr<Vector<Real>> &control,
82  const Ptr<Vector<Real>> &adjoint,
83  const Ptr<Vector<Real>> &dualstate,
84  const Ptr<Vector<Real>> &dualcontrol,
85  const Ptr<Vector<Real>> &dualadjoint,
86  const bool storage,
87  const bool useFDhessVec)
88  : obj_(obj), con_(con),
89  storage_(storage), useFDhessVec_(useFDhessVec),
90  nupda_(0), nvalu_(0), ngrad_(0), nhess_(0), nprec_(0),
91  nstat_(0), nadjo_(0), nssen_(0), nasen_(0),
92  updateFlag_(true), updateIter_(0), updateType_(UpdateType::Initial),
93  newUpdate_(false) {
94  stateStore_ = makePtr<VectorController<Real>>();
95  adjointStore_ = makePtr<VectorController<Real>>();
96  state_ = state->clone(); state_->set(*state);
97  adjoint_ = adjoint->clone();
98  state_sens_ = state->clone();
99  adjoint_sens_ = adjoint->clone();
100  dualstate_ = dualstate->clone();
101  dualstate1_ = dualstate->clone();
102  dualadjoint_ = dualadjoint->clone();
103  dualcontrol_ = dualcontrol->clone();
104 }
105 
106 template<typename Real>
108  const Ptr<Objective_SimOpt<Real>> &obj,
109  const Ptr<Constraint_SimOpt<Real>> &con,
110  const Ptr<VectorController<Real>> &stateStore,
111  const Ptr<Vector<Real>> &state,
112  const Ptr<Vector<Real>> &control,
113  const Ptr<Vector<Real>> &adjoint,
114  const bool storage,
115  const bool useFDhessVec)
116  : obj_(obj), con_(con), stateStore_(stateStore),
117  storage_(storage), useFDhessVec_(useFDhessVec),
118  nupda_(0), nvalu_(0), ngrad_(0), nhess_(0), nprec_(0),
119  nstat_(0), nadjo_(0), nssen_(0), nasen_(0),
120  updateFlag_(true), updateIter_(0), updateType_(UpdateType::Initial),
121  newUpdate_(false) {
122  adjointStore_ = makePtr<VectorController<Real>>();
123  state_ = state->clone(); state_->set(*state);
124  adjoint_ = adjoint->clone();
125  state_sens_ = state->clone();
126  adjoint_sens_ = adjoint->clone();
127  dualstate_ = state->dual().clone();
128  dualstate1_ = state->dual().clone();
129  dualadjoint_ = adjoint->dual().clone();
130  dualcontrol_ = control->dual().clone();
131 }
132 
133 template<typename Real>
135  const Ptr<Objective_SimOpt<Real>> &obj,
136  const Ptr<Constraint_SimOpt<Real>> &con,
137  const Ptr<VectorController<Real>> &stateStore,
138  const Ptr<Vector<Real>> &state,
139  const Ptr<Vector<Real>> &control,
140  const Ptr<Vector<Real>> &adjoint,
141  const Ptr<Vector<Real>> &dualstate,
142  const Ptr<Vector<Real>> &dualcontrol,
143  const Ptr<Vector<Real>> &dualadjoint,
144  const bool storage,
145  const bool useFDhessVec)
146  : obj_(obj), con_(con), stateStore_(stateStore),
147  storage_(storage), useFDhessVec_(useFDhessVec),
148  nupda_(0), nvalu_(0), ngrad_(0), nhess_(0), nprec_(0),
149  nstat_(0), nadjo_(0), nssen_(0), nasen_(0),
150  updateFlag_(true), updateIter_(0), updateType_(UpdateType::Initial),
151  newUpdate_(false) {
152  adjointStore_ = makePtr<VectorController<Real>>();
153  state_ = state->clone(); state_->set(*state);
154  adjoint_ = adjoint->clone();
155  state_sens_ = state->clone();
156  adjoint_sens_ = adjoint->clone();
157  dualstate_ = dualstate->clone();
158  dualstate1_ = dualstate->clone();
159  dualadjoint_ = dualadjoint->clone();
160  dualcontrol_ = dualcontrol->clone();
161 }
162 
163 template<typename Real>
164 void Reduced_Objective_SimOpt<Real>::update( const Vector<Real> &z, bool flag, int iter ) {
165  nupda_++;
166  newUpdate_ = false;
167  updateFlag_ = flag;
168  updateIter_ = iter;
169  stateStore_->objectiveUpdate(true);
170  adjointStore_->objectiveUpdate(flag);
171 }
172 
173 template<typename Real>
175  nupda_++;
176  newUpdate_ = true;
177  updateType_ = type;
178  updateIter_ = iter;
179  stateStore_->objectiveUpdate(type);
180  adjointStore_->objectiveUpdate(type);
181 }
182 
183 template<typename Real>
185  nvalu_++;
186  // Solve state equation
187  solve_state_equation(z,tol);
188  // Get objective function value
189  return obj_->value(*state_,z,tol);
190 }
191 
192 template<typename Real>
194  ngrad_++;
195  // Solve state equation
196  solve_state_equation(z,tol);
197  // Solve adjoint equation
198  solve_adjoint_equation(z,tol);
199  // Evaluate the full gradient wrt z
200  obj_->gradient_2(*dualcontrol_,*state_,z,tol);
201  // Build gradient
202  con_->applyAdjointJacobian_2(g,*adjoint_,*state_,z,tol);
203  g.plus(*dualcontrol_);
204 }
205 
206 template<typename Real>
208  nhess_++;
209  if ( useFDhessVec_ ) {
210  Objective<Real>::hessVec(hv,v,z,tol);
211  }
212  else {
213  // Solve state equation
214  solve_state_equation(z,tol);
215  // Solve adjoint equation
216  solve_adjoint_equation(z,tol);
217  // Solve state sensitivity equation
218  solve_state_sensitivity(v,z,tol);
219  // Solve adjoint sensitivity equation
220  solve_adjoint_sensitivity(v,z,tol);
221  // Build hessVec
222  con_->applyAdjointJacobian_2(hv,*adjoint_sens_,*state_,z,tol);
223  obj_->hessVec_21(*dualcontrol_,*state_sens_,*state_,z,tol);
224  hv.plus(*dualcontrol_);
225  obj_->hessVec_22(*dualcontrol_,v,*state_,z,tol);
226  hv.plus(*dualcontrol_);
227  con_->applyAdjointHessian_12(*dualcontrol_,*adjoint_,*state_sens_,*state_,z,tol);
228  hv.plus(*dualcontrol_);
229  con_->applyAdjointHessian_22(*dualcontrol_,*adjoint_,v,*state_,z,tol);
230  hv.plus(*dualcontrol_);
231  }
232 }
233 
234 template<typename Real>
236  nprec_++;
237  Pv.set(v.dual());
238 }
239 
240 template<typename Real>
241 void Reduced_Objective_SimOpt<Real>::setParameter(const std::vector<Real> &param) {
243  con_->setParameter(param);
244  obj_->setParameter(param);
245 }
246 
247 template<typename Real>
248 void Reduced_Objective_SimOpt<Real>::summarize(std::ostream &stream) const {
249  stream << std::endl;
250  stream << std::string(80,'=') << std::endl;
251  stream << " ROL::Reduced_Objective_SimOpt::summarize" << std::endl;
252  stream << " Number of calls to update: " << nupda_ << std::endl;
253  stream << " Number of calls to value: " << nvalu_ << std::endl;
254  stream << " Number of calls to gradient: " << ngrad_ << std::endl;
255  stream << " Number of calls to hessvec: " << nhess_ << std::endl;
256  stream << " Number of calls to precond: " << nprec_ << std::endl;
257  stream << " Number of state solves: " << nstat_ << std::endl;
258  stream << " Number of adjoint solves: " << nadjo_ << std::endl;
259  stream << " Number of state sensitivity solves: " << nssen_ << std::endl;
260  stream << " Number of adjoint sensitivity solves: " << nssen_ << std::endl;
261  stream << std::string(80,'=') << std::endl;
262  stream << std::endl;
263 }
264 
265 template<typename Real>
267  nupda_ = 0; nvalu_ = 0; ngrad_ = 0; nhess_ = 0; nprec_ = 0;
268  nstat_ = 0; nadjo_ = 0; nssen_ = 0; nasen_ = 0;
269 }
270 
271 template<typename Real>
273  // Check if state has been computed.
274  bool isComputed = storage_ ? stateStore_->get(*state_,Objective<Real>::getParameter()) : false;
275  // Solve state equation if not done already.
276  if (!isComputed || !storage_) {
277  // Update equality constraint with new Opt variable.
278  if (newUpdate_) con_->update_2(z,updateType_,updateIter_);
279  else con_->update_2(z,updateFlag_,updateIter_);
280  // Solve state equation.
281  con_->solve(*dualadjoint_,*state_,z,tol);
282  nstat_++;
283  // Update equality constraint with new Sim variable.
284  if (newUpdate_) con_->update_1(*state_,updateType_,updateIter_);
285  else con_->update_1(*state_,updateFlag_,updateIter_);
286  // Update full objective function.
287  if (newUpdate_) obj_->update(*state_,z,updateType_,updateIter_);
288  else obj_->update(*state_,z,updateFlag_,updateIter_);
289  // Store state.
290  if (storage_) stateStore_->set(*state_,Objective<Real>::getParameter());
291  }
292 }
293 
294 template<typename Real>
296  // Check if adjoint has been computed.
297  bool isComputed = storage_ ? adjointStore_->get(*adjoint_,Objective<Real>::getParameter()) : false;
298  // Solve adjoint equation if not done already.
299  if (!isComputed || !storage_) {
300  // Evaluate the full gradient wrt u
301  obj_->gradient_1(*dualstate_,*state_,z,tol);
302  // Solve adjoint equation
303  con_->applyInverseAdjointJacobian_1(*adjoint_,*dualstate_,*state_,z,tol);
304  adjoint_->scale(static_cast<Real>(-1));
305  nadjo_++;
306  // Store adjoint
307  if (storage_) adjointStore_->set(*adjoint_,Objective<Real>::getParameter());
308  }
309 }
310 
311 template<typename Real>
313  // Solve state sensitivity equation
314  con_->applyJacobian_2(*dualadjoint_,v,*state_,z,tol);
315  dualadjoint_->scale(static_cast<Real>(-1));
316  con_->applyInverseJacobian_1(*state_sens_,*dualadjoint_,*state_,z,tol);
317  nssen_++;
318 }
319 
320 template<typename Real>
322  // Evaluate full hessVec in the direction (s,v)
323  obj_->hessVec_11(*dualstate_,*state_sens_,*state_,z,tol);
324  obj_->hessVec_12(*dualstate1_,v,*state_,z,tol);
325  dualstate_->plus(*dualstate1_);
326  // Apply adjoint Hessian of constraint
327  con_->applyAdjointHessian_11(*dualstate1_,*adjoint_,*state_sens_,*state_,z,tol);
328  dualstate_->plus(*dualstate1_);
329  con_->applyAdjointHessian_21(*dualstate1_,*adjoint_,v,*state_,z,tol);
330  dualstate_->plus(*dualstate1_);
331  // Solve adjoint sensitivity equation
332  dualstate_->scale(static_cast<Real>(-1));
333  con_->applyInverseAdjointJacobian_1(*adjoint_sens_,*dualstate_,*state_,z,tol);
334  nasen_++;
335 }
336 
337 } // namespace ROL
338 
339 #endif
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 hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
void solve_state_equation(const Vector< Real > &z, Real &tol)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
Ptr< VectorController< Real > > adjointStore_
virtual void precond(Vector< Real > &Pv, const Vector< Real > &v, const Vector< Real > &z, Real &tol) override
Apply a reduced Hessian preconditioner.
void solve_adjoint_sensitivity(const Vector< Real > &v, const Vector< Real > &z, Real &tol)
Given , the adjoint variable , and a direction , solve the adjoint sensitvity equation for ...
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:226
Ptr< VectorController< Real > > stateStore_
Real value(const Vector< Real > &z, Real &tol) override
Given , evaluate the objective function where solves .
void update(const Vector< Real > &z, bool flag=true, int iter=-1) override
Update the SimOpt objective function and equality constraint.
void solve_adjoint_equation(const Vector< Real > &z, Real &tol)
Given which solves the state equation, solve the adjoint equation for .
void summarize(std::ostream &stream) const
virtual void setParameter(const std::vector< Real > &param)
void setParameter(const std::vector< Real > &param) override
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
void gradient(Vector< Real > &g, const Vector< Real > &z, Real &tol) override
Given , evaluate the gradient of the objective function where solves .
Defines the constraint operator interface for simulation-based optimization.
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &z, Real &tol) override
Given , evaluate the Hessian of the objective function in the direction .
Reduced_Objective_SimOpt(const Ptr< Objective_SimOpt< Real >> &obj, const Ptr< Constraint_SimOpt< Real >> &con, const Ptr< Vector< Real >> &state, const Ptr< Vector< Real >> &control, const Ptr< Vector< Real >> &adjoint, const bool storage=true, const bool useFDhessVec=false)
Constructor.
void solve_state_sensitivity(const Vector< Real > &v, const Vector< Real > &z, Real &tol)
Given which solves the state equation and a direction , solve the state senstivity equation for ...
const Ptr< Obj > obj_