ROL
ROL_CylinderHead.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
49#ifndef USE_HESSVEC
50#define USE_HESSVEC 1
51#endif
52
53#ifndef ROL_CYLINDERHEAD_HPP
54#define ROL_CYLINDERHEAD_HPP
55
57#include "ROL_StdObjective.hpp"
58#include "ROL_StdConstraint.hpp"
59#include "ROL_TestProblem.hpp"
60
61namespace ROL {
62namespace ZOO {
63
64 template<class Real>
65 class Objective_CylinderHead : public StdObjective<Real> {
66 private:
67 Real warranty(const Real flatness, const int deriv = 0) const {
68 const Real w0(1e5), w1(1.5e4), w2(4), zero(0);
69 return (deriv==0 ? w0 + w1*(w2 - flatness)
70 :(deriv==1 ? -w1 : zero));
71 }
72
73 Real horsepower(const Real dintake, const int deriv = 0) const {
74 const Real d0(250), d1(200), d2(1), d3(1.833), zero(0);
75 return (deriv==0 ? d0 + d1*(dintake/d3 - d2)
76 :(deriv==1 ? d1/d3 : zero));
77 }
78
79 public:
81
82 Real value( const std::vector<Real> &x, Real &tol ) {
83 const Real w0(1e5), h0(250);
84 return -(horsepower(x[0],0)/h0 + warranty(x[1],0)/w0);
85 }
86
87 void gradient( std::vector<Real> &g, const std::vector<Real> &x, Real &tol ) {
88 const Real w0(1e5), h0(250);
89 g[0] = -horsepower(x[0],1)/h0;
90 g[1] = -warranty(x[1],1)/w0;
91 }
92#if USE_HESSVEC
93 void hessVec( std::vector<Real> &hv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol ) {
94 const Real w0(1e5), h0(250);
95 hv[0] = -horsepower(x[0],2)/h0 * v[0];
96 hv[1] = -warranty(x[1],2)/w0 * v[1];
97 }
98#endif
99 };
100
101 template<class Real>
103 private:
104 Real warranty(const Real flatness, const int deriv = 0) const {
105 const Real w0(1e5), w1(1.5e4), w2(4), zero(0);
106 return (deriv==0 ? w0 + w1*(w2 - flatness)
107 :(deriv==1 ? -w1 : zero));
108 }
109
110 Real tcycle(const Real flatness, const int deriv = 0) const {
111 const Real t0(45), t1(4.5), t2(4), tpwr(1.5), one(1);
112 return (deriv==0 ? t0 + t1*std::pow(t2 - flatness, tpwr)
113 :(deriv==1 ? -t1*tpwr*std::sqrt(t2 - flatness)
114 : t1*tpwr*(tpwr-one)/std::sqrt(t2 - flatness)));
115 }
116
117 Real twall(const Real dintake, const int deriv = 0) const {
118 const Real ointake(3.25), oexhaust(1.34), dexhaust(1.556), half(0.5), zero(0);
119 return (deriv==0 ? ointake - oexhaust - half*(dintake + dexhaust)
120 :(deriv==1 ? -half : zero));
121 }
122
123 Real smax(const Real tw, const int deriv = 0) const {
124 const Real s0(750), spwr(2.5), one(1), two(2);
125 return (deriv==0 ? s0 + one/std::pow(tw, spwr)
126 :(deriv==1 ? -spwr/std::pow(tw, spwr+one)
127 : spwr*(spwr+one)/std::pow(tw, spwr+two)));
128 }
129
130 public:
132
133 void value( std::vector<Real> &c, const std::vector<Real> &x, Real &tol ) {
134 const Real one(1), two(2), sixty(60), syield(3000), w0(1e5);
135 c[0] = two*smax(twall(x[0],0),0)/syield - one;
136 c[1] = one - warranty(x[1],0)/w0;
137 c[2] = tcycle(x[1],0)/sixty - one;
138 }
139
140 void applyJacobian( std::vector<Real> &jv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol ) {
141 const Real two(2), sixty(60), syield(3000), w0(1e5);
142 jv[0] = two*smax(twall(x[0],0),1)/syield * twall(x[0],1) * v[0];
143 jv[1] = -warranty(x[1],1)/w0 * v[1];
144 jv[2] = tcycle(x[1],1)/sixty * v[1];
145 }
146
147 void applyAdjointJacobian( std::vector<Real> &ajv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol ) {
148 const Real two(2), sixty(60), syield(3000), w0(1e5);
149 ajv[0] = two*smax(twall(x[0],0),1)/syield * twall(x[0],1) * v[0];
150 ajv[1] = -warranty(x[1],1)/w0 * v[1] + tcycle(x[1],1)/sixty * v[2];
151 }
152#if USE_HESSVEC
153 void applyAdjointHessian( std::vector<Real> &ahuv, const std::vector<Real> &u, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol ) {
154 const Real two(2), sixty(60), syield(3000), w0(1e5);
155 Real tw = twall(x[0],0);
156 ahuv[0] = two*(smax(tw,2) * twall(x[0],1) * twall(x[0],1)
157 + smax(tw,1) * twall(x[0],2))/syield * u[0] * v[0];
158 ahuv[1] = (-warranty(x[1],2)/w0 * u[1] + tcycle(x[1],2)/sixty * u[2]) * v[1];
159 }
160#endif
161 };
162
163 template<class Real>
164 class getCylinderHead : public TestProblem<Real> {
165 public:
167
168 Ptr<Objective<Real>> getObjective(void) const {
169 return makePtr<Objective_CylinderHead<Real>>();
170 }
171
172 Ptr<Vector<Real>> getInitialGuess(void) const {
173 int n = 2;
174 Ptr<std::vector<Real>> scale = makePtr<std::vector<Real>>(n,static_cast<Real>(1.0));
175 Ptr<std::vector<Real>> xp = makePtr<std::vector<Real>>(n,static_cast<Real>(0.0));
176 (*xp)[0] = static_cast<Real>(1.8);
177 (*xp)[1] = static_cast<Real>(1.0);
178 return makePtr<PrimalScaledStdVector<Real>>(xp,scale);
179 }
180
181 Ptr<Vector<Real>> getSolution(const int i = 0) const {
182 int n = 2;
183 Ptr<std::vector<Real>> scale = makePtr<std::vector<Real>>(n,static_cast<Real>(1.0));
184 Ptr<std::vector<Real>> xp = makePtr<std::vector<Real>>(n,static_cast<Real>(0.0));
185 (*xp)[0] = static_cast<Real>(2.122);
186 (*xp)[1] = static_cast<Real>(1.769);
187 return makePtr<PrimalScaledStdVector<Real>>(xp,scale);
188 }
189
190 Ptr<BoundConstraint<Real>> getBoundConstraint(void) const {
191 int n = 2;
192 Ptr<std::vector<Real>> scale = makePtr<std::vector<Real>>(n,static_cast<Real>(1.0));
193 Ptr<std::vector<Real>> lp = makePtr<std::vector<Real>>(n,static_cast<Real>(0.0));
194 Ptr<std::vector<Real>> up = makePtr<std::vector<Real>>(n,static_cast<Real>(0.0));
195 (*lp)[0] = static_cast<Real>(1.5);
196 (*lp)[1] = static_cast<Real>(0.0);
197 (*up)[0] = static_cast<Real>(2.164);
198 (*up)[1] = static_cast<Real>(4.0);
199 Ptr<Vector<Real>> l = makePtr<PrimalScaledStdVector<Real>>(lp,scale);
200 Ptr<Vector<Real>> u = makePtr<PrimalScaledStdVector<Real>>(up,scale);
201 return makePtr<Bounds<Real>>(l,u);
202 }
203
204 Ptr<Constraint<Real>> getInequalityConstraint(void) const {
205 return makePtr<Constraint_CylinderHead<Real>>();
206 }
207
208 Ptr<Vector<Real>> getInequalityMultiplier(void) const {
209 Ptr<std::vector<Real>> scale = makePtr<std::vector<Real>>(3,static_cast<Real>(1.0));
210 Ptr<std::vector<Real>> lp = makePtr<std::vector<Real>>(3,static_cast<Real>(0.0));
211 return makePtr<DualScaledStdVector<Real>>(lp,scale);
212 }
213
214 Ptr<BoundConstraint<Real>> getSlackBoundConstraint(void) const {
215 Ptr<std::vector<Real>> scale = makePtr<std::vector<Real>>(3,static_cast<Real>(1.0));
216 Ptr<std::vector<Real>> up = makePtr<std::vector<Real>>(3,static_cast<Real>(0.0));
217 Ptr<Vector<Real>> u = makePtr<DualScaledStdVector<Real>>(up,scale);
218 return makePtr<Bounds<Real>>(*u,false);
219 }
220 };
221
222}// End ZOO Namespace
223}// End ROL Namespace
224
225#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Contains definitions of test objective functions.
Defines the equality constraint operator interface for StdVectors.
void applyAdjointHessian(Vector< Real > &ahuv, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &x, Real &tol) override
Apply the derivative of the adjoint of the constraint Jacobian at to vector in direction ,...
Specializes the ROL::Objective interface for objective functions that operate on ROL::StdVector's.
virtual void hessVec(std::vector< Real > &hv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
Real tcycle(const Real flatness, const int deriv=0) const
void value(std::vector< Real > &c, const std::vector< Real > &x, Real &tol)
Real warranty(const Real flatness, const int deriv=0) const
Real twall(const Real dintake, const int deriv=0) const
void applyAdjointJacobian(std::vector< Real > &ajv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
void applyJacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
Real smax(const Real tw, const int deriv=0) const
void gradient(std::vector< Real > &g, const std::vector< Real > &x, Real &tol)
Real warranty(const Real flatness, const int deriv=0) const
Real value(const std::vector< Real > &x, Real &tol)
Real horsepower(const Real dintake, const int deriv=0) const
Ptr< Vector< Real > > getInitialGuess(void) const
Ptr< Vector< Real > > getInequalityMultiplier(void) const
Ptr< BoundConstraint< Real > > getBoundConstraint(void) const
Ptr< Vector< Real > > getSolution(const int i=0) const
Ptr< Constraint< Real > > getInequalityConstraint(void) const
Ptr< BoundConstraint< Real > > getSlackBoundConstraint(void) const
Ptr< Objective< Real > > getObjective(void) const