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
47namespace ROL {
48
49template<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), isUpdated_(true) {
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
76template<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), isUpdated_(true) {
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
106template<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), isUpdated_(true) {
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
133template<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), isUpdated_(true) {
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
163template<typename Real>
164void Reduced_Objective_SimOpt<Real>::update( const Vector<Real> &z, bool flag, int iter ) {
165 nupda_++;
166 isUpdated_ = false;
167 newUpdate_ = false;
168 updateFlag_ = flag;
169 updateIter_ = iter;
170 stateStore_->objectiveUpdate(true);
171 adjointStore_->objectiveUpdate(flag);
172}
173
174template<typename Real>
176 nupda_++;
177 isUpdated_ = false;
178 newUpdate_ = true;
179 updateType_ = type;
180 updateIter_ = iter;
181 stateStore_->objectiveUpdate(type);
182 adjointStore_->objectiveUpdate(type);
183}
184
185template<typename Real>
187 nvalu_++;
188 // Solve state equation
189 solve_state_equation(z,tol);
190 // Get objective function value
191 return obj_->value(*state_,z,tol);
192}
193
194template<typename Real>
196 ngrad_++;
197 // Solve state equation
198 solve_state_equation(z,tol);
199 // Solve adjoint equation
200 solve_adjoint_equation(z,tol);
201 // Evaluate the full gradient wrt z
202 obj_->gradient_2(*dualcontrol_,*state_,z,tol);
203 // Build gradient
204 con_->applyAdjointJacobian_2(g,*adjoint_,*state_,z,tol);
205 g.plus(*dualcontrol_);
206}
207
208template<typename Real>
210 nhess_++;
211 if ( useFDhessVec_ ) {
212 Objective<Real>::hessVec(hv,v,z,tol);
213 }
214 else {
215 // Solve state equation
216 solve_state_equation(z,tol);
217 // Solve adjoint equation
218 solve_adjoint_equation(z,tol);
219 // Solve state sensitivity equation
220 solve_state_sensitivity(v,z,tol);
221 // Solve adjoint sensitivity equation
222 solve_adjoint_sensitivity(v,z,tol);
223 // Build hessVec
224 con_->applyAdjointJacobian_2(hv,*adjoint_sens_,*state_,z,tol);
225 obj_->hessVec_21(*dualcontrol_,*state_sens_,*state_,z,tol);
226 hv.plus(*dualcontrol_);
227 obj_->hessVec_22(*dualcontrol_,v,*state_,z,tol);
228 hv.plus(*dualcontrol_);
229 con_->applyAdjointHessian_12(*dualcontrol_,*adjoint_,*state_sens_,*state_,z,tol);
230 hv.plus(*dualcontrol_);
231 con_->applyAdjointHessian_22(*dualcontrol_,*adjoint_,v,*state_,z,tol);
232 hv.plus(*dualcontrol_);
233 }
234}
235
236template<typename Real>
238 nprec_++;
239 Pv.set(v.dual());
240}
241
242template<typename Real>
243void Reduced_Objective_SimOpt<Real>::setParameter(const std::vector<Real> &param) {
245 con_->setParameter(param);
246 obj_->setParameter(param);
247}
248
249template<typename Real>
250void Reduced_Objective_SimOpt<Real>::summarize(std::ostream &stream, const Ptr<BatchManager<Real>> &bman) const {
251 int nupda(0), nvalu(0), ngrad(0), nhess(0), nprec(0), nstat(0), nadjo(0), nssen(0), nasen(0);
252 if (bman == nullPtr) {
253 nupda = nupda_;
254 nvalu = nvalu_;
255 ngrad = ngrad_;
256 nhess = nhess_;
257 nprec = nprec_;
258 nstat = nstat_;
259 nadjo = nadjo_;
260 nssen = nssen_;
261 nasen = nasen_;
262 }
263 else {
264 auto sumAll = [bman](int val) {
265 Real global(0), local(val);
266 bman->sumAll(&local,&global,1);
267 return static_cast<int>(global);
268 };
269 nupda = sumAll(nupda_);
270 nvalu = sumAll(nvalu_);
271 ngrad = sumAll(ngrad_);
272 nhess = sumAll(nhess_);
273 nprec = sumAll(nprec_);
274 nstat = sumAll(nstat_);
275 nadjo = sumAll(nadjo_);
276 nssen = sumAll(nssen_);
277 nasen = sumAll(nasen_);
278 }
279 stream << std::endl;
280 stream << std::string(80,'=') << std::endl;
281 stream << " ROL::Reduced_Objective_SimOpt::summarize" << std::endl;
282 stream << " Number of calls to update: " << nupda << std::endl;
283 stream << " Number of calls to value: " << nvalu << std::endl;
284 stream << " Number of calls to gradient: " << ngrad << std::endl;
285 stream << " Number of calls to hessvec: " << nhess << std::endl;
286 stream << " Number of calls to precond: " << nprec << std::endl;
287 stream << " Number of state solves: " << nstat << std::endl;
288 stream << " Number of adjoint solves: " << nadjo << std::endl;
289 stream << " Number of state sensitivity solves: " << nssen << std::endl;
290 stream << " Number of adjoint sensitivity solves: " << nasen << std::endl;
291 stream << std::string(80,'=') << std::endl;
292 stream << std::endl;
293}
294
295template<typename Real>
297 nupda_ = 0; nvalu_ = 0; ngrad_ = 0; nhess_ = 0; nprec_ = 0;
298 nstat_ = 0; nadjo_ = 0; nssen_ = 0; nasen_ = 0;
299}
300
301template<typename Real>
303 if (!isUpdated_) {
304 // Update equality constraint with new Opt variable.
305 if (newUpdate_) con_->update_2(z,updateType_,updateIter_);
306 else con_->update_2(z,updateFlag_,updateIter_);
307 }
308 // Check if state has been computed.
309 bool isComputed = storage_ ? stateStore_->get(*state_,Objective<Real>::getParameter()) : false;
310 // Solve state equation if not done already.
311 if (!isComputed || !storage_) {
312 // Solve state equation.
313 con_->solve(*dualadjoint_,*state_,z,tol);
314 nstat_++;
315 // Store state.
316 if (storage_) stateStore_->set(*state_,Objective<Real>::getParameter());
317 }
318 if (!isUpdated_) {
319 // Update equality constraint with new Sim variable.
320 if (newUpdate_) con_->update_1(*state_,updateType_,updateIter_);
321 else con_->update_1(*state_,updateFlag_,updateIter_);
322 // Update full objective function.
323 if (newUpdate_) obj_->update(*state_,z,updateType_,updateIter_);
324 else obj_->update(*state_,z,updateFlag_,updateIter_);
325 isUpdated_ = true;
326 }
327}
328
329template<typename Real>
331 // Check if adjoint has been computed.
332 bool isComputed = storage_ ? adjointStore_->get(*adjoint_,Objective<Real>::getParameter()) : false;
333 // Solve adjoint equation if not done already.
334 if (!isComputed || !storage_) {
335 // Evaluate the full gradient wrt u
336 obj_->gradient_1(*dualstate_,*state_,z,tol);
337 // Solve adjoint equation
338 con_->applyInverseAdjointJacobian_1(*adjoint_,*dualstate_,*state_,z,tol);
339 adjoint_->scale(static_cast<Real>(-1));
340 nadjo_++;
341 // Store adjoint
342 if (storage_) adjointStore_->set(*adjoint_,Objective<Real>::getParameter());
343 }
344}
345
346template<typename Real>
348 // Solve state sensitivity equation
349 con_->applyJacobian_2(*dualadjoint_,v,*state_,z,tol);
350 dualadjoint_->scale(static_cast<Real>(-1));
351 con_->applyInverseJacobian_1(*state_sens_,*dualadjoint_,*state_,z,tol);
352 nssen_++;
353}
354
355template<typename Real>
357 // Evaluate full hessVec in the direction (s,v)
358 obj_->hessVec_11(*dualstate_,*state_sens_,*state_,z,tol);
359 obj_->hessVec_12(*dualstate1_,v,*state_,z,tol);
360 dualstate_->plus(*dualstate1_);
361 // Apply adjoint Hessian of constraint
362 con_->applyAdjointHessian_11(*dualstate1_,*adjoint_,*state_sens_,*state_,z,tol);
363 dualstate_->plus(*dualstate1_);
364 con_->applyAdjointHessian_21(*dualstate1_,*adjoint_,v,*state_,z,tol);
365 dualstate_->plus(*dualstate1_);
366 // Solve adjoint sensitivity equation
367 dualstate_->scale(static_cast<Real>(-1));
368 con_->applyInverseAdjointJacobian_1(*adjoint_sens_,*dualstate_,*state_,z,tol);
369 nasen_++;
370}
371
372} // namespace ROL
373
374#endif
const Ptr< Obj > obj_
Defines the constraint operator interface for simulation-based optimization.
Provides the interface to evaluate simulation-based objective functions.
Provides the interface to evaluate objective functions.
virtual 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.
virtual void precond(Vector< Real > &Pv, const Vector< Real > &v, const Vector< Real > &z, Real &tol) override
Apply a reduced Hessian preconditioner.
Ptr< VectorController< Real > > adjointStore_
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 .
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 .
Real value(const Vector< Real > &z, Real &tol) override
Given , evaluate the objective function where solves .
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_adjoint_equation(const Vector< Real > &z, Real &tol)
Given which solves the state equation, solve the adjoint equation for .
Ptr< VectorController< Real > > stateStore_
void update(const Vector< Real > &z, bool flag=true, int iter=-1) override
Update the SimOpt objective function and equality constraint.
void solve_state_equation(const Vector< Real > &z, Real &tol)
void gradient(Vector< Real > &g, const Vector< Real > &z, Real &tol) override
Given , evaluate the gradient of the objective function where solves .
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 .
void setParameter(const std::vector< Real > &param) override
void summarize(std::ostream &stream, const Ptr< BatchManager< Real > > &bman=nullPtr) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
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
virtual void plus(const Vector &x)=0
Compute , where .