ROL
ROL_SerialConstraint.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#pragma once
45#ifndef ROL_SERIALCONSTRAINT_HPP
46#define ROL_SERIALCONSTRAINT_HPP
47
48#include <type_traits>
49
53
61namespace ROL {
62
63template<typename Real>
65 public SerialFunction<Real> {
66private:
67
69 using SerialFunction<Real>::ts;
70 using SerialFunction<Real>::clone;
71
72 Ptr<DynamicConstraint<Real>> con_; // Constraint over a single time step
73
74public:
75
76 using size_type = typename std::vector<Real>::size_type;
81
83 const Vector<Real>& u_initial,
84 const TimeStampsPtr<Real>& timeStampsPtr ) :
85 SerialFunction<Real>::SerialFunction( u_initial, timeStampsPtr ),
86 con_(con) {}
87
88 virtual void solve( Vector<Real>& c,
89 Vector<Real>& u,
90 const Vector<Real>& z,
91 Real& tol ) override {
92
93 auto& cp = partition(c);
94 auto& up = partition(u);
95 auto& zp = partition(z);
96
98 con_->solve( cp[0], getInitialCondition(), up[0], zp[0], ts(0) );
99
100 for( size_type k=1; k<numTimeSteps(); ++k )
101 con_->solve( cp[k], up[k-1], up[k], zp[k], ts(k) );
102
103 } // solve
104
105 virtual void update( const Vector<Real> &u, const Vector<Real> &z, bool flag = true, int iter = -1 ) override {
106 auto& up = partition(u);
107 auto& zp = partition(z);
108
110 con_->update( getInitialCondition(), up[0], zp[0], ts(0) );
111
112 for( size_type k=1; k<numTimeSteps(); ++k )
113 con_->update( up[k-1], up[k], zp[k], ts(k) );
114 }
115
116
117 using Constraint_SimOpt<Real>::value;
118 virtual void value( Vector<Real>& c,
119 const Vector<Real>& u,
120 const Vector<Real>& z,
121 Real& tol ) override {
122
123 auto& cp = partition(c);
124 auto& up = partition(u);
125 auto& zp = partition(z);
126
128 con_->value( cp[0], getInitialCondition(), up[0], zp[0], ts(0) );
129
130 for( size_type k=1; k<numTimeSteps(); ++k )
131 con_->value( cp[k], up[k-1], up[k], zp[k], ts(k) );
132
133 } // value
134
135
136 virtual void applyJacobian_1( Vector<Real>& jv,
137 const Vector<Real>& v,
138 const Vector<Real>& u,
139 const Vector<Real>& z,
140 Real& tol ) override {
141
142 auto& jvp = partition(jv); auto& vp = partition(v);
143 auto& up = partition(u); auto& zp = partition(z);
144
145 auto tmp = clone(jvp[0]);
146 auto& x = *tmp;
147
149 con_->applyJacobian_un( jvp[0], vp[0], getZeroState(), up[0], zp[0], ts(0) );
150
151 for( size_type k=1; k<numTimeSteps(); ++k ) {
152
153 con_->applyJacobian_uo( x, vp[k-1], up[k-1], up[k], zp[k], ts(k) );
154 con_->applyJacobian_un( jvp[k], vp[k], up[k-1], up[k], zp[k], ts(k) );
155 jvp.get(k)->plus(x);
156
157 } // end for
158
159 } // applyJacobian_1
160
161
163 const Vector<Real>& v,
164 const Vector<Real>& u,
165 const Vector<Real>& z,
166 Real& tol) override {
167
168 auto& ijvp = partition(ijv); auto& vp = partition(v);
169 auto& up = partition(u); auto& zp = partition(z);
170
171 auto tmp = clone(ijvp[0]);
172 auto& x = *tmp;
173
175 con_->applyInverseJacobian_un( ijvp[0], vp[0], getZeroState(),
176 up[0], zp[0], ts(0) );
177
178 for( size_type k=1; k<numTimeSteps(); ++k ) {
179
180 con_->applyJacobian_uo( x, ijvp[k-1], up[k-1], up[k], zp[k], ts(k) );
181 x.scale(-1.0);
182 x.plus( vp[k] );
183 con_->applyInverseJacobian_un( ijvp[k], x, up[k-1], up[k], zp[k], ts(k) );
184
185 } // end for
186
187 } // applyInverseJacobian_1
188
189
192 const Vector<Real>& v,
193 const Vector<Real>& u,
194 const Vector<Real>& z,
195 const Vector<Real>& dualv,
196 Real& tol) override {
197
198 auto& ajvp = partition(ajv); auto& vp = partition(v);
199 auto& up = partition(u); auto& zp = partition(z);
200
201 auto tmp = clone(ajvp[0]);
202 auto& x = *tmp;
203
205 con_->applyAdjointJacobian_un( ajvp[0], vp[0], getZeroState(), up[0], zp[0], ts(0) );
206
207 for( size_type k=1; k<numTimeSteps(); ++k ) {
208
209 con_->applyAdjointJacobian_un( ajvp[k], vp[k], up[k-1], up[k], zp[k], ts(k) );
210 con_->applyAdjointJacobian_uo( x, vp[k], up[k-1], up[k], zp[k], ts(k) );
211 ajvp[k-1].plus(x);
212
213 } // end for
214
215 } // applyAdjointJacobian_1
216
218 const Vector<Real>& v,
219 const Vector<Real>& u,
220 const Vector<Real>& z,
221 Real& tol) override {
222
223 auto& iajvp = partition(iajv); auto& vp = partition(v);
224 auto& up = partition(u); auto& zp = partition(z);
225
226 auto tmp = clone(iajvp[0]);
227 auto& x = *tmp;
228
229 size_type k = numTimeSteps()-1;
230
231 con_->applyInverseAdjointJacobian_un( iajvp[k], vp[k], up[k-1], up[k], zp[k], ts(k) );
232
233 for( size_type k=numTimeSteps()-2; k>0; --k ) {
234
235 con_->applyAdjointJacobian_uo( x, iajvp[k+1], up[k], up[k+1], zp[k+1], ts(k+1) );
236 x.scale(-1.0);
237 x.plus( vp[k] );
238 con_->applyInverseAdjointJacobian_un( iajvp[k], x, up[k-1], up[k], zp[k], ts(k) );
239
240 } // end for
241
242 con_->applyAdjointJacobian_uo( x, iajvp[1], up[0], up[1], zp[1], ts(1) );
243 x.scale(-1.0);
244 x.plus( vp[0] );
245
247 con_->applyInverseAdjointJacobian_un( iajvp[0], x, getZeroState(), up[0], zp[0], ts(0) );
248
249 // this weird condition places iajvp in the final vector slot
250 else iajvp[0].set(x);
251
252 } // applyInverseAdjointJacobian_1
253
254
255 virtual void applyJacobian_2( Vector<Real>& jv,
256 const Vector<Real>& v,
257 const Vector<Real>& u,
258 const Vector<Real>& z,
259 Real &tol ) override {
260
261 auto& jvp = partition(jv); auto& vp = partition(v);
262 auto& up = partition(u); auto& zp = partition(z);
263
265 con_->applyJacobian_z( jvp[0], vp[0], getInitialCondition(), up[0], zp[0], ts(0) );
266
267 for( size_type k=1; k<numTimeSteps(); ++k )
268 con_->applyJacobian_z( jvp[k], vp[k], up[k-1], up[k], zp[k], ts(k) );
269
270 } // applyJacobian_2
271
272
275 const Vector<Real>& v,
276 const Vector<Real>& u,
277 const Vector<Real>& z,
278 Real& tol ) override {
279
280 auto& ajvp = partition(ajv); auto& vp = partition(v);
281 auto& up = partition(u); auto& zp = partition(z);
282
284 con_->applyAdjointJacobian_z( ajvp[0], vp[0], getInitialCondition(), up[0], zp[0], ts(0) );
285
286 for( size_type k=1; k<numTimeSteps(); ++k )
287 con_->applyAdjointJacobian_z( ajvp[k], vp[k], up[k-1], up[k], zp[k], ts(k) );
288
289 } // applyAdjointJacobian_2
290
291
292
294 const Vector<Real>& w,
295 const Vector<Real>& v,
296 const Vector<Real>& u,
297 const Vector<Real>& z,
298 Real& tol) override {
299
300 auto& ahwvp = partition(ahwv); auto& wp = partition(w);
301 auto& vp = partition(v); auto& up = partition(u);
302 auto& zp = partition(z);
303
304 auto tmp = clone(ahwvp[0]);
305 auto& x = *tmp;
306
307 if( !getSkipInitialCondition() ) {
308
309 con_->applyAdjointHessian_un_un( ahwvp[0], wp[0], vp[0], getZeroState(), up[0], zp[0], ts(0) );
310 con_->applyAdjointHessian_un_uo( x, wp[1], vp[1], up[0], up[1], zp[1], ts(1) );
311 ahwvp[0].plus(x);
312 con_->applyAdjointHessian_uo_uo( x, wp[1], vp[0], up[0], up[1], zp[1], ts(1) );
313 ahwvp[0].plus(x);
314
315 }
316
317 for( size_type k=1; k<numTimeSteps(); ++k ) {
318
319 con_->applyAdjointHessian_un_un( ahwvp[k], wp[k], vp[k], up[k-1], up[k], zp[k], ts(k) );
320 con_->applyAdjointHessian_uo_un( x, wp[k], vp[k-1], up[k-1], up[k], zp[k], ts(k) );
321 ahwvp[k].plus(x);
322
323 if( k < numTimeSteps()-1 ) {
324 con_->applyAdjointHessian_un_uo( x, wp[k+1], vp[k+1], up[k], up[k+1], zp[k+1], ts(k+1) );
325 ahwvp[k].plus(x);
326 con_->applyAdjointHessian_uo_uo( x, wp[k+1], vp[k], up[k], up[k+1], zp[k+1], ts(k+1) );
327 ahwvp[k].plus(x);
328 } // endif
329 } // endfor
330
331 } // applyAdjointHessian_11
332
333
334 // TODO: Implement off-diagonal blocks
336 const Vector<Real>& w,
337 const Vector<Real>& v,
338 const Vector<Real>& u,
339 const Vector<Real>& z,
340 Real& tol) override {
341 }
342
343
345 const Vector<Real>& w,
346 const Vector<Real>& v,
347 const Vector<Real>& u,
348 const Vector<Real>& z,
349 Real& tol) override {
350 }
351
353 const Vector<Real>& w,
354 const Vector<Real>& v,
355 const Vector<Real>& u,
356 const Vector<Real>& z,
357 Real& tol) override {
358
359 auto& ahwvp = partition(ahwv); auto& vp = partition(v); auto& wp = partition(w);
360 auto& up = partition(u); auto& zp = partition(z);
361
362 con_->applyAdjointHessian_z_z( ahwvp[0], wp[0], vp[0], getInitialCondition(),
363 up[0], zp[0], ts(0) );
364
365 for( size_type k=1; k<numTimeSteps(); ++k ) {
366 con_->applyAdjointHessian_z_z( ahwvp[k], wp[k], vp[k], up[k-1],
367 up[k], zp[k], ts(k) );
368 }
369
370 } // applyAdjointHessian_22
371
372
373}; // SerialConstraint
374
375// Helper function to create a new SerialConstraint
376
377template<typename DynCon, typename Real, typename P = Ptr<SerialConstraint<Real>>>
378inline typename std::enable_if<std::is_base_of<DynamicConstraint<Real>,DynCon>::value,P>::type
379make_SerialConstraint( const Ptr<DynCon>& con,
380 const Vector<Real>& u_initial,
381 const TimeStampsPtr<Real> timeStampsPtr ) {
382 return makePtr<SerialConstraint<Real>>(con,u_initial,timeStampsPtr);
383}
384
385
386
387
388
389
390} // namespace ROL
391
392
393#endif // ROL_SERIALCONSTRAINT_HPP
Defines the constraint operator interface for simulation-based optimization.
Defines the time-dependent constraint operator interface for simulation-based optimization.
Defines the linear algebra of vector space on a generic partitioned vector.
Evaluates ROL::DynamicConstraint over a sequential set of time intervals.
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...
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 ...
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 .
SerialConstraint(const Ptr< DynamicConstraint< Real > > &con, const Vector< Real > &u_initial, const TimeStampsPtr< Real > &timeStampsPtr)
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 .
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 .
typename std::vector< Real >::size_type size_type
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 solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Given , solve for .
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 i...
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 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 ...
Ptr< DynamicConstraint< Real > > con_
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...
virtual void value(Vector< Real > &c, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Evaluate the constraint operator at .
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...
Provides behavior common to SerialObjective as SerialConstaint.
size_type numTimeSteps() const
const Vector< Real > & getInitialCondition() const
const Vector< Real > & getZeroState() const
bool getSkipInitialCondition() const
Ptr< Vector< Real > > clone(const Vector< Real > &x)
const TimeStamp< Real > & ts(size_type i) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
Ptr< std::vector< TimeStamp< Real > > > TimeStampsPtr
ROL::Objective_SerialSimOpt Objective_SimOpt value(const V &u, const V &z, Real &tol) override
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)
PartitionedVector< Real > & partition(Vector< Real > &x)