ROL
ROL_SerialConstraint.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 #pragma once
11 #ifndef ROL_SERIALCONSTRAINT_HPP
12 #define ROL_SERIALCONSTRAINT_HPP
13 
14 #include <type_traits>
15 
18 #include "ROL_SerialFunction.hpp"
19 
27 namespace ROL {
28 
29 template<typename Real>
30 class SerialConstraint : public Constraint_SimOpt<Real>,
31  public SerialFunction<Real> {
32 private:
33 
37 
38  Ptr<DynamicConstraint<Real>> con_; // Constraint over a single time step
39 
40 public:
41 
47 
49  const Vector<Real>& u_initial,
50  const TimeStampsPtr<Real>& timeStampsPtr ) :
51  SerialFunction<Real>::SerialFunction( u_initial, timeStampsPtr ),
52  con_(con) {}
53 
54  virtual void solve( Vector<Real>& c,
55  Vector<Real>& u,
56  const Vector<Real>& z,
57  Real& tol ) override {
58 
59  auto& cp = partition(c);
60  auto& up = partition(u);
61  auto& zp = partition(z);
62 
64  con_->solve( cp[0], getInitialCondition(), up[0], zp[0], ts(0) );
65 
66  for( size_type k=1; k<numTimeSteps(); ++k )
67  con_->solve( cp[k], up[k-1], up[k], zp[k], ts(k) );
68 
69  } // solve
70 
71  virtual void update( const Vector<Real> &u, const Vector<Real> &z, bool flag = true, int iter = -1 ) override {
72  auto& up = partition(u);
73  auto& zp = partition(z);
74 
76  con_->update( getInitialCondition(), up[0], zp[0], ts(0) );
77 
78  for( size_type k=1; k<numTimeSteps(); ++k )
79  con_->update( up[k-1], up[k], zp[k], ts(k) );
80  }
81 
82 
84  virtual void value( Vector<Real>& c,
85  const Vector<Real>& u,
86  const Vector<Real>& z,
87  Real& tol ) override {
88 
89  auto& cp = partition(c);
90  auto& up = partition(u);
91  auto& zp = partition(z);
92 
94  con_->value( cp[0], getInitialCondition(), up[0], zp[0], ts(0) );
95 
96  for( size_type k=1; k<numTimeSteps(); ++k )
97  con_->value( cp[k], up[k-1], up[k], zp[k], ts(k) );
98 
99  } // value
100 
101 
102  virtual void applyJacobian_1( Vector<Real>& jv,
103  const Vector<Real>& v,
104  const Vector<Real>& u,
105  const Vector<Real>& z,
106  Real& tol ) override {
107 
108  auto& jvp = partition(jv); auto& vp = partition(v);
109  auto& up = partition(u); auto& zp = partition(z);
110 
111  auto tmp = clone(jvp[0]);
112  auto& x = *tmp;
113 
114  if( !getSkipInitialCondition() )
115  con_->applyJacobian_un( jvp[0], vp[0], getZeroState(), up[0], zp[0], ts(0) );
116 
117  for( size_type k=1; k<numTimeSteps(); ++k ) {
118 
119  con_->applyJacobian_uo( x, vp[k-1], up[k-1], up[k], zp[k], ts(k) );
120  con_->applyJacobian_un( jvp[k], vp[k], up[k-1], up[k], zp[k], ts(k) );
121  jvp.get(k)->plus(x);
122 
123  } // end for
124 
125  } // applyJacobian_1
126 
127 
129  const Vector<Real>& v,
130  const Vector<Real>& u,
131  const Vector<Real>& z,
132  Real& tol) override {
133 
134  auto& ijvp = partition(ijv); auto& vp = partition(v);
135  auto& up = partition(u); auto& zp = partition(z);
136 
137  auto tmp = clone(ijvp[0]);
138  auto& x = *tmp;
139 
140  if( !getSkipInitialCondition() )
141  con_->applyInverseJacobian_un( ijvp[0], vp[0], getZeroState(),
142  up[0], zp[0], ts(0) );
143 
144  for( size_type k=1; k<numTimeSteps(); ++k ) {
145 
146  con_->applyJacobian_uo( x, ijvp[k-1], up[k-1], up[k], zp[k], ts(k) );
147  x.scale(-1.0);
148  x.plus( vp[k] );
149  con_->applyInverseJacobian_un( ijvp[k], x, up[k-1], up[k], zp[k], ts(k) );
150 
151  } // end for
152 
153  } // applyInverseJacobian_1
154 
155 
158  const Vector<Real>& v,
159  const Vector<Real>& u,
160  const Vector<Real>& z,
161  const Vector<Real>& dualv,
162  Real& tol) override {
163 
164  auto& ajvp = partition(ajv); auto& vp = partition(v);
165  auto& up = partition(u); auto& zp = partition(z);
166 
167  auto tmp = clone(ajvp[0]);
168  auto& x = *tmp;
169 
170  if( !getSkipInitialCondition() )
171  con_->applyAdjointJacobian_un( ajvp[0], vp[0], getZeroState(), up[0], zp[0], ts(0) );
172 
173  for( size_type k=1; k<numTimeSteps(); ++k ) {
174 
175  con_->applyAdjointJacobian_un( ajvp[k], vp[k], up[k-1], up[k], zp[k], ts(k) );
176  con_->applyAdjointJacobian_uo( x, vp[k], up[k-1], up[k], zp[k], ts(k) );
177  ajvp[k-1].plus(x);
178 
179  } // end for
180 
181  } // applyAdjointJacobian_1
182 
184  const Vector<Real>& v,
185  const Vector<Real>& u,
186  const Vector<Real>& z,
187  Real& tol) override {
188 
189  auto& iajvp = partition(iajv); auto& vp = partition(v);
190  auto& up = partition(u); auto& zp = partition(z);
191 
192  auto tmp = clone(iajvp[0]);
193  auto& x = *tmp;
194 
195  size_type k = numTimeSteps()-1;
196 
197  con_->applyInverseAdjointJacobian_un( iajvp[k], vp[k], up[k-1], up[k], zp[k], ts(k) );
198 
199  for( size_type k=numTimeSteps()-2; k>0; --k ) {
200 
201  con_->applyAdjointJacobian_uo( x, iajvp[k+1], up[k], up[k+1], zp[k+1], ts(k+1) );
202  x.scale(-1.0);
203  x.plus( vp[k] );
204  con_->applyInverseAdjointJacobian_un( iajvp[k], x, up[k-1], up[k], zp[k], ts(k) );
205 
206  } // end for
207 
208  con_->applyAdjointJacobian_uo( x, iajvp[1], up[0], up[1], zp[1], ts(1) );
209  x.scale(-1.0);
210  x.plus( vp[0] );
211 
212  if( !getSkipInitialCondition() )
213  con_->applyInverseAdjointJacobian_un( iajvp[0], x, getZeroState(), up[0], zp[0], ts(0) );
214 
215  // this weird condition places iajvp in the final vector slot
216  else iajvp[0].set(x);
217 
218  } // applyInverseAdjointJacobian_1
219 
220 
221  virtual void applyJacobian_2( Vector<Real>& jv,
222  const Vector<Real>& v,
223  const Vector<Real>& u,
224  const Vector<Real>& z,
225  Real &tol ) override {
226 
227  auto& jvp = partition(jv); auto& vp = partition(v);
228  auto& up = partition(u); auto& zp = partition(z);
229 
230  if( !getSkipInitialCondition() )
231  con_->applyJacobian_z( jvp[0], vp[0], getInitialCondition(), up[0], zp[0], ts(0) );
232 
233  for( size_type k=1; k<numTimeSteps(); ++k )
234  con_->applyJacobian_z( jvp[k], vp[k], up[k-1], up[k], zp[k], ts(k) );
235 
236  } // applyJacobian_2
237 
238 
241  const Vector<Real>& v,
242  const Vector<Real>& u,
243  const Vector<Real>& z,
244  Real& tol ) override {
245 
246  auto& ajvp = partition(ajv); auto& vp = partition(v);
247  auto& up = partition(u); auto& zp = partition(z);
248 
249  if( !getSkipInitialCondition() )
250  con_->applyAdjointJacobian_z( ajvp[0], vp[0], getInitialCondition(), up[0], zp[0], ts(0) );
251 
252  for( size_type k=1; k<numTimeSteps(); ++k )
253  con_->applyAdjointJacobian_z( ajvp[k], vp[k], up[k-1], up[k], zp[k], ts(k) );
254 
255  } // applyAdjointJacobian_2
256 
257 
258 
259  virtual void applyAdjointHessian_11( Vector<Real>& ahwv,
260  const Vector<Real>& w,
261  const Vector<Real>& v,
262  const Vector<Real>& u,
263  const Vector<Real>& z,
264  Real& tol) override {
265 
266  auto& ahwvp = partition(ahwv); auto& wp = partition(w);
267  auto& vp = partition(v); auto& up = partition(u);
268  auto& zp = partition(z);
269 
270  auto tmp = clone(ahwvp[0]);
271  auto& x = *tmp;
272 
273  if( !getSkipInitialCondition() ) {
274 
275  con_->applyAdjointHessian_un_un( ahwvp[0], wp[0], vp[0], getZeroState(), up[0], zp[0], ts(0) );
276  con_->applyAdjointHessian_un_uo( x, wp[1], vp[1], up[0], up[1], zp[1], ts(1) );
277  ahwvp[0].plus(x);
278  con_->applyAdjointHessian_uo_uo( x, wp[1], vp[0], up[0], up[1], zp[1], ts(1) );
279  ahwvp[0].plus(x);
280 
281  }
282 
283  for( size_type k=1; k<numTimeSteps(); ++k ) {
284 
285  con_->applyAdjointHessian_un_un( ahwvp[k], wp[k], vp[k], up[k-1], up[k], zp[k], ts(k) );
286  con_->applyAdjointHessian_uo_un( x, wp[k], vp[k-1], up[k-1], up[k], zp[k], ts(k) );
287  ahwvp[k].plus(x);
288 
289  if( k < numTimeSteps()-1 ) {
290  con_->applyAdjointHessian_un_uo( x, wp[k+1], vp[k+1], up[k], up[k+1], zp[k+1], ts(k+1) );
291  ahwvp[k].plus(x);
292  con_->applyAdjointHessian_uo_uo( x, wp[k+1], vp[k], up[k], up[k+1], zp[k+1], ts(k+1) );
293  ahwvp[k].plus(x);
294  } // endif
295  } // endfor
296 
297  } // applyAdjointHessian_11
298 
299 
300  // TODO: Implement off-diagonal blocks
302  const Vector<Real>& w,
303  const Vector<Real>& v,
304  const Vector<Real>& u,
305  const Vector<Real>& z,
306  Real& tol) override {
307  }
308 
309 
311  const Vector<Real>& w,
312  const Vector<Real>& v,
313  const Vector<Real>& u,
314  const Vector<Real>& z,
315  Real& tol) override {
316  }
317 
318  virtual void applyAdjointHessian_22( Vector<Real>& ahwv,
319  const Vector<Real>& w,
320  const Vector<Real>& v,
321  const Vector<Real>& u,
322  const Vector<Real>& z,
323  Real& tol) override {
324 
325  auto& ahwvp = partition(ahwv); auto& vp = partition(v); auto& wp = partition(w);
326  auto& up = partition(u); auto& zp = partition(z);
327 
328  con_->applyAdjointHessian_z_z( ahwvp[0], wp[0], vp[0], getInitialCondition(),
329  up[0], zp[0], ts(0) );
330 
331  for( size_type k=1; k<numTimeSteps(); ++k ) {
332  con_->applyAdjointHessian_z_z( ahwvp[k], wp[k], vp[k], up[k-1],
333  up[k], zp[k], ts(k) );
334  }
335 
336  } // applyAdjointHessian_22
337 
338 
339 }; // SerialConstraint
340 
341 // Helper function to create a new SerialConstraint
342 
343 template<typename DynCon, typename Real, typename P = Ptr<SerialConstraint<Real>>>
344 inline typename std::enable_if<std::is_base_of<DynamicConstraint<Real>,DynCon>::value,P>::type
345 make_SerialConstraint( const Ptr<DynCon>& con,
346  const Vector<Real>& u_initial,
347  const TimeStampsPtr<Real> timeStampsPtr ) {
348  return makePtr<SerialConstraint<Real>>(con,u_initial,timeStampsPtr);
349 }
350 
351 
352 
353 
354 
355 
356 } // namespace ROL
357 
358 
359 #endif // ROL_SERIALCONSTRAINT_HPP
PartitionedVector< Real > & partition(Vector< Real > &x)
const Vector< Real > & getInitialCondition() const
typename PV< Real >::size_type size_type
virtual void applyAdjointHessian_12(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the optimization-space derivative of the adjoint of the constraint simulation-space Jacobian at...
const Vector< Real > & getZeroState() const
Defines the time-dependent constraint operator interface for simulation-based optimization.
Defines the linear algebra of vector space on a generic partitioned vector.
virtual void applyAdjointHessian_21(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the simulation-space derivative of the adjoint of the constraint optimization-space Jacobian at...
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1) override
Update constraint functions. x is the optimization variable, flag = true if optimization variable is ...
virtual void applyAdjointHessian_22(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the optimization-space derivative of the adjoint of the constraint optimization-space Jacobian ...
ROL::Objective_SimOpt value
size_type numTimeSteps() const
virtual void value(Vector< Real > &c, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Evaluate the constraint operator at .
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual void applyAdjointJacobian_1(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualv, Real &tol) override
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the secondary int...
virtual void applyJacobian_1(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the partial constraint Jacobian at , , to the vector .
virtual void applyJacobian_2(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the partial constraint Jacobian at , , to the vector .
std::enable_if< std::is_base_of< DynamicConstraint< Real >, DynCon >::value, P >::type make_SerialConstraint(const Ptr< DynCon > &con, const Vector< Real > &u_initial, const TimeStampsPtr< Real > timeStampsPtr)
Ptr< Vector< Real > > clone(const Vector< Real > &x)
typename std::vector< Real >::size_type size_type
virtual void applyAdjointJacobian_2(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the primary interface...
SerialConstraint(const Ptr< DynamicConstraint< Real >> &con, const Vector< Real > &u_initial, const TimeStampsPtr< Real > &timeStampsPtr)
void applyInverseAdjointJacobian_1(Vector< Real > &iajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the inverse of the adjoint of the partial constraint Jacobian at , , to the vector ...
bool getSkipInitialCondition() const
Provides behavior common to SerialObjective as SerialConstaint.
virtual void applyAdjointHessian_11(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the simulation-space derivative of the adjoint of the constraint simulation-space Jacobian at ...
Ptr< DynamicConstraint< Real > > con_
Ptr< std::vector< TimeStamp< Real >>> TimeStampsPtr
virtual void solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Given , solve for .
Evaluates ROL::DynamicConstraint over a sequential set of time intervals.
Defines the constraint operator interface for simulation-based optimization.
const TimeStamp< Real > & ts(size_type i) const
virtual void applyInverseJacobian_1(Vector< Real > &ijv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the inverse partial constraint Jacobian at , , to the vector .