ROL
ROL_Reduced_Objective_SimOpt_Def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_REDUCED_OBJECTIVE_SIMOPT_DEF_H
11 #define ROL_REDUCED_OBJECTIVE_SIMOPT_DEF_H
12 
13 namespace ROL {
14 
15 template<typename Real>
17  const Ptr<Objective_SimOpt<Real>> &obj,
18  const Ptr<Constraint_SimOpt<Real>> &con,
19  const Ptr<Vector<Real>> &state,
20  const Ptr<Vector<Real>> &control,
21  const Ptr<Vector<Real>> &adjoint,
22  const bool storage,
23  const bool useFDhessVec)
24  : obj_(obj), con_(con),
25  storage_(storage), useFDhessVec_(useFDhessVec),
26  nupda_(0), nvalu_(0), ngrad_(0), nhess_(0), nprec_(0),
27  nstat_(0), nadjo_(0), nssen_(0), nasen_(0),
28  updateFlag_(true), updateIter_(0), updateType_(UpdateType::Initial),
29  newUpdate_(false), isUpdated_(true) {
30  stateStore_ = makePtr<VectorController<Real>>();
31  adjointStore_ = makePtr<VectorController<Real>>();
32  state_ = state->clone(); state_->set(*state);
33  adjoint_ = adjoint->clone();
34  state_sens_ = state->clone();
35  adjoint_sens_ = adjoint->clone();
36  dualstate_ = state->dual().clone();
37  dualstate1_ = state->dual().clone();
38  dualadjoint_ = adjoint->dual().clone();
39  dualcontrol_ = control->dual().clone();
40 }
41 
42 template<typename Real>
44  const Ptr<Objective_SimOpt<Real>> &obj,
45  const Ptr<Constraint_SimOpt<Real>> &con,
46  const Ptr<Vector<Real>> &state,
47  const Ptr<Vector<Real>> &control,
48  const Ptr<Vector<Real>> &adjoint,
49  const Ptr<Vector<Real>> &dualstate,
50  const Ptr<Vector<Real>> &dualcontrol,
51  const Ptr<Vector<Real>> &dualadjoint,
52  const bool storage,
53  const bool useFDhessVec)
54  : obj_(obj), con_(con),
55  storage_(storage), useFDhessVec_(useFDhessVec),
56  nupda_(0), nvalu_(0), ngrad_(0), nhess_(0), nprec_(0),
57  nstat_(0), nadjo_(0), nssen_(0), nasen_(0),
58  updateFlag_(true), updateIter_(0), updateType_(UpdateType::Initial),
59  newUpdate_(false), isUpdated_(true) {
60  stateStore_ = makePtr<VectorController<Real>>();
61  adjointStore_ = makePtr<VectorController<Real>>();
62  state_ = state->clone(); state_->set(*state);
63  adjoint_ = adjoint->clone();
64  state_sens_ = state->clone();
65  adjoint_sens_ = adjoint->clone();
66  dualstate_ = dualstate->clone();
67  dualstate1_ = dualstate->clone();
68  dualadjoint_ = dualadjoint->clone();
69  dualcontrol_ = dualcontrol->clone();
70 }
71 
72 template<typename Real>
74  const Ptr<Objective_SimOpt<Real>> &obj,
75  const Ptr<Constraint_SimOpt<Real>> &con,
76  const Ptr<VectorController<Real>> &stateStore,
77  const Ptr<Vector<Real>> &state,
78  const Ptr<Vector<Real>> &control,
79  const Ptr<Vector<Real>> &adjoint,
80  const bool storage,
81  const bool useFDhessVec)
82  : obj_(obj), con_(con), stateStore_(stateStore),
83  storage_(storage), useFDhessVec_(useFDhessVec),
84  nupda_(0), nvalu_(0), ngrad_(0), nhess_(0), nprec_(0),
85  nstat_(0), nadjo_(0), nssen_(0), nasen_(0),
86  updateFlag_(true), updateIter_(0), updateType_(UpdateType::Initial),
87  newUpdate_(false), isUpdated_(true) {
88  adjointStore_ = makePtr<VectorController<Real>>();
89  state_ = state->clone(); state_->set(*state);
90  adjoint_ = adjoint->clone();
91  state_sens_ = state->clone();
92  adjoint_sens_ = adjoint->clone();
93  dualstate_ = state->dual().clone();
94  dualstate1_ = state->dual().clone();
95  dualadjoint_ = adjoint->dual().clone();
96  dualcontrol_ = control->dual().clone();
97 }
98 
99 template<typename Real>
101  const Ptr<Objective_SimOpt<Real>> &obj,
102  const Ptr<Constraint_SimOpt<Real>> &con,
103  const Ptr<VectorController<Real>> &stateStore,
104  const Ptr<Vector<Real>> &state,
105  const Ptr<Vector<Real>> &control,
106  const Ptr<Vector<Real>> &adjoint,
107  const Ptr<Vector<Real>> &dualstate,
108  const Ptr<Vector<Real>> &dualcontrol,
109  const Ptr<Vector<Real>> &dualadjoint,
110  const bool storage,
111  const bool useFDhessVec)
112  : obj_(obj), con_(con), stateStore_(stateStore),
113  storage_(storage), useFDhessVec_(useFDhessVec),
114  nupda_(0), nvalu_(0), ngrad_(0), nhess_(0), nprec_(0),
115  nstat_(0), nadjo_(0), nssen_(0), nasen_(0),
116  updateFlag_(true), updateIter_(0), updateType_(UpdateType::Initial),
117  newUpdate_(false), isUpdated_(true) {
118  adjointStore_ = makePtr<VectorController<Real>>();
119  state_ = state->clone(); state_->set(*state);
120  adjoint_ = adjoint->clone();
121  state_sens_ = state->clone();
122  adjoint_sens_ = adjoint->clone();
123  dualstate_ = dualstate->clone();
124  dualstate1_ = dualstate->clone();
125  dualadjoint_ = dualadjoint->clone();
126  dualcontrol_ = dualcontrol->clone();
127 }
128 
129 template<typename Real>
130 void Reduced_Objective_SimOpt<Real>::update( const Vector<Real> &z, bool flag, int iter ) {
131  nupda_++;
132  isUpdated_ = false;
133  newUpdate_ = false;
134  updateFlag_ = flag;
135  updateIter_ = iter;
136  stateStore_->objectiveUpdate(true);
137  adjointStore_->objectiveUpdate(flag);
138 }
139 
140 template<typename Real>
142  nupda_++;
143  isUpdated_ = false;
144  newUpdate_ = true;
145  updateType_ = type;
146  updateIter_ = iter;
147  stateStore_->objectiveUpdate(type);
148  adjointStore_->objectiveUpdate(type);
149 }
150 
151 template<typename Real>
153  nvalu_++;
154  // Solve state equation
155  solve_state_equation(z,tol);
156  // Get objective function value
157  return obj_->value(*state_,z,tol);
158 }
159 
160 template<typename Real>
162  ngrad_++;
163  // Solve state equation
164  solve_state_equation(z,tol);
165  // Solve adjoint equation
166  solve_adjoint_equation(z,tol);
167  // Evaluate the full gradient wrt z
168  obj_->gradient_2(*dualcontrol_,*state_,z,tol);
169  // Build gradient
170  con_->applyAdjointJacobian_2(g,*adjoint_,*state_,z,tol);
171  g.plus(*dualcontrol_);
172 }
173 
174 template<typename Real>
176  nhess_++;
177  if ( useFDhessVec_ ) {
178  Objective<Real>::hessVec(hv,v,z,tol);
179  }
180  else {
181  // Solve state equation
182  solve_state_equation(z,tol);
183  // Solve adjoint equation
184  solve_adjoint_equation(z,tol);
185  // Solve state sensitivity equation
186  solve_state_sensitivity(v,z,tol);
187  // Solve adjoint sensitivity equation
188  solve_adjoint_sensitivity(v,z,tol);
189  // Build hessVec
190  con_->applyAdjointJacobian_2(hv,*adjoint_sens_,*state_,z,tol);
191  obj_->hessVec_21(*dualcontrol_,*state_sens_,*state_,z,tol);
192  hv.plus(*dualcontrol_);
193  obj_->hessVec_22(*dualcontrol_,v,*state_,z,tol);
194  hv.plus(*dualcontrol_);
195  con_->applyAdjointHessian_12(*dualcontrol_,*adjoint_,*state_sens_,*state_,z,tol);
196  hv.plus(*dualcontrol_);
197  con_->applyAdjointHessian_22(*dualcontrol_,*adjoint_,v,*state_,z,tol);
198  hv.plus(*dualcontrol_);
199  }
200 }
201 
202 template<typename Real>
204  nprec_++;
205  Pv.set(v.dual());
206 }
207 
208 template<typename Real>
209 void Reduced_Objective_SimOpt<Real>::setParameter(const std::vector<Real> &param) {
211  con_->setParameter(param);
212  obj_->setParameter(param);
213 }
214 
215 template<typename Real>
216 void Reduced_Objective_SimOpt<Real>::summarize(std::ostream &stream, const Ptr<BatchManager<Real>> &bman) const {
217  int nupda(0), nvalu(0), ngrad(0), nhess(0), nprec(0), nstat(0), nadjo(0), nssen(0), nasen(0);
218  if (bman == nullPtr) {
219  nupda = nupda_;
220  nvalu = nvalu_;
221  ngrad = ngrad_;
222  nhess = nhess_;
223  nprec = nprec_;
224  nstat = nstat_;
225  nadjo = nadjo_;
226  nssen = nssen_;
227  nasen = nasen_;
228  }
229  else {
230  auto sumAll = [bman](int val) {
231  Real global(0), local(val);
232  bman->sumAll(&local,&global,1);
233  return static_cast<int>(global);
234  };
235  nupda = sumAll(nupda_);
236  nvalu = sumAll(nvalu_);
237  ngrad = sumAll(ngrad_);
238  nhess = sumAll(nhess_);
239  nprec = sumAll(nprec_);
240  nstat = sumAll(nstat_);
241  nadjo = sumAll(nadjo_);
242  nssen = sumAll(nssen_);
243  nasen = sumAll(nasen_);
244  }
245  stream << std::endl;
246  stream << std::string(80,'=') << std::endl;
247  stream << " ROL::Reduced_Objective_SimOpt::summarize" << std::endl;
248  stream << " Number of calls to update: " << nupda << std::endl;
249  stream << " Number of calls to value: " << nvalu << std::endl;
250  stream << " Number of calls to gradient: " << ngrad << std::endl;
251  stream << " Number of calls to hessvec: " << nhess << std::endl;
252  stream << " Number of calls to precond: " << nprec << std::endl;
253  stream << " Number of state solves: " << nstat << std::endl;
254  stream << " Number of adjoint solves: " << nadjo << std::endl;
255  stream << " Number of state sensitivity solves: " << nssen << std::endl;
256  stream << " Number of adjoint sensitivity solves: " << nasen << std::endl;
257  stream << std::string(80,'=') << std::endl;
258  stream << std::endl;
259 }
260 
261 template<typename Real>
263  nupda_ = 0; nvalu_ = 0; ngrad_ = 0; nhess_ = 0; nprec_ = 0;
264  nstat_ = 0; nadjo_ = 0; nssen_ = 0; nasen_ = 0;
265 }
266 
267 template<typename Real>
269  if (!isUpdated_) {
270  // Update equality constraint with new Opt variable.
271  if (newUpdate_) con_->update_2(z,updateType_,updateIter_);
272  else con_->update_2(z,updateFlag_,updateIter_);
273  }
274  // Check if state has been computed.
275  bool isComputed = storage_ ? stateStore_->get(*state_,Objective<Real>::getParameter()) : false;
276  // Solve state equation if not done already.
277  if (!isComputed || !storage_) {
278  // Solve state equation.
279  con_->solve(*dualadjoint_,*state_,z,tol);
280  nstat_++;
281  // Store state.
282  if (storage_) stateStore_->set(*state_,Objective<Real>::getParameter());
283  }
284  if (!isUpdated_) {
285  // Update equality constraint with new Sim variable.
286  if (newUpdate_) con_->update_1(*state_,updateType_,updateIter_);
287  else con_->update_1(*state_,updateFlag_,updateIter_);
288  // Update full objective function.
289  if (newUpdate_) obj_->update(*state_,z,updateType_,updateIter_);
290  else obj_->update(*state_,z,updateFlag_,updateIter_);
291  isUpdated_ = true;
292  }
293 }
294 
295 template<typename Real>
297  // Check if adjoint has been computed.
298  bool isComputed = storage_ ? adjointStore_->get(*adjoint_,Objective<Real>::getParameter()) : false;
299  // Solve adjoint equation if not done already.
300  if (!isComputed || !storage_) {
301  // Evaluate the full gradient wrt u
302  obj_->gradient_1(*dualstate_,*state_,z,tol);
303  // Solve adjoint equation
304  con_->applyInverseAdjointJacobian_1(*adjoint_,*dualstate_,*state_,z,tol);
305  adjoint_->scale(static_cast<Real>(-1));
306  nadjo_++;
307  // Store adjoint
308  if (storage_) adjointStore_->set(*adjoint_,Objective<Real>::getParameter());
309  }
310 }
311 
312 template<typename Real>
314  // Solve state sensitivity equation
315  con_->applyJacobian_2(*dualadjoint_,v,*state_,z,tol);
316  dualadjoint_->scale(static_cast<Real>(-1));
317  con_->applyInverseJacobian_1(*state_sens_,*dualadjoint_,*state_,z,tol);
318  nssen_++;
319 }
320 
321 template<typename Real>
323  // Evaluate full hessVec in the direction (s,v)
324  obj_->hessVec_11(*dualstate_,*state_sens_,*state_,z,tol);
325  obj_->hessVec_12(*dualstate1_,v,*state_,z,tol);
326  dualstate_->plus(*dualstate1_);
327  // Apply adjoint Hessian of constraint
328  con_->applyAdjointHessian_11(*dualstate1_,*adjoint_,*state_sens_,*state_,z,tol);
329  dualstate_->plus(*dualstate1_);
330  con_->applyAdjointHessian_21(*dualstate1_,*adjoint_,v,*state_,z,tol);
331  dualstate_->plus(*dualstate1_);
332  // Solve adjoint sensitivity equation
333  dualstate_->scale(static_cast<Real>(-1));
334  con_->applyInverseAdjointJacobian_1(*adjoint_sens_,*dualstate_,*state_,z,tol);
335  nasen_++;
336 }
337 
338 } // namespace ROL
339 
340 #endif
Provides the interface to evaluate objective functions.
Provides the interface to evaluate simulation-based objective functions.
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:192
virtual void plus(const Vector &x)=0
Compute , where .
void summarize(std::ostream &stream, const Ptr< BatchManager< Real >> &bman=nullPtr) const
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:46
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 ...
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 .
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:175
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_