1 #pragma once
2
3 #include "ROL_StdVector.hpp"
4 #include "ROL_Objective_SimOpt.hpp"
5 #include "ROL_Constraint_SimOpt.hpp"
6
7 namespace Rocket {
8
getVector(ROL::Vector<double> & x)9 std::vector<double>& getVector( ROL::Vector<double>& x ) {
10 return *(dynamic_cast<ROL::StdVector<double>&>(x).getVector());
11 }
12
getVector(const ROL::Vector<double> & x)13 const std::vector<double>& getVector( const ROL::Vector<double>& x ) {
14 return *(dynamic_cast<const ROL::StdVector<double>&>(x).getVector());
15 }
16
17
18 class Objective : public ROL::Objective_SimOpt<double> {
19 private:
20 using V = ROL::Vector<double>;
21
22 int N;
23 double T, dt, mt;
24 double htarg, alpha;
25 const ROL::Ptr<const V> w; // Trapezoidal weights
26
27 public:
28
Objective(int N_,double T_,double mt_,double htarg_,double alpha_,const ROL::Ptr<const V> & w_)29 Objective( int N_, double T_, double mt_, double htarg_, double alpha_,
30 const ROL::Ptr<const V>& w_ ) :
31 N(N_), T(T_), dt(T/N), mt(mt_), htarg(htarg_), alpha(alpha_), w(w_) {
32 }
33
value(const V & u,const V & z,double & tol)34 double value( const V& u, const V& z, double& tol ) {
35 return 0.5*std::pow(htarg-w->dot(u),2) + 0.5*alpha*dt*z.dot(z);
36 }
37
gradient_1(V & g,const V & u,const V & z,double & tol)38 void gradient_1( V& g, const V& u, const V& z, double& tol ) {
39 g.set(*w); g.scale(w->dot(u)-htarg);
40 }
41
gradient_2(V & g,const V & u,const V & z,double & tol)42 void gradient_2( V& g, const V& u, const V& z, double& tol ) {
43 g.set(z); g.scale(alpha*dt);
44 }
45
hessVec_11(V & hv,const V & v,const V & u,const V & z,double & tol)46 void hessVec_11( V& hv, const V& v, const V& u, const V& z, double& tol ) {
47 hv.set(*w); hv.scale(w->dot(v));
48 }
49
hessVec_12(V & hv,const V & v,const V & u,const V & z,double & tol)50 void hessVec_12( V& hv, const V& v, const V& u, const V& z, double& tol ) {
51 hv.zero();
52 }
53
hessVec_21(V & hv,const V & v,const V & u,const V & z,double & tol)54 void hessVec_21( V& hv, const V& v, const V& u, const V& z, double& tol ) {
55 hv.zero();
56 }
57
hessVec_22(V & hv,const V & v,const V & u,const V & z,double & tol)58 void hessVec_22( V& hv, const V& v, const V& u, const V& z, double& tol ) {
59 hv.set(v); hv.scale(alpha*dt);
60 }
61 }; // Objective
62
63
64
65
66 class Constraint : public ROL::Constraint_SimOpt<double> {
67 private:
68 using V = ROL::Vector<double>;
69
70 int N;
71 double T, dt, gdt, mt, mf, ve;
72
73 std::vector<double> mass;
74
75 public:
76
Constraint(int N_,double T_,double mt_,double mf_,double ve_,double g_)77 Constraint( int N_, double T_, double mt_,
78 double mf_, double ve_, double g_ ) : N(N_), T(T_),
79 dt(T/N), gdt(g_*dt), mt(mt_),
80 mf(mf_), ve(ve_), mass(N) {
81 mass[0] = mt;
82 }
83
update_2(const V & z,bool flag=true,int iter=-1)84 void update_2( const V& z, bool flag = true, int iter = -1 ) {
85 auto& zs = getVector(z);
86
87 mass[0] = mt - dt*zs[0];
88 for( int i=1; i<N; ++i )
89 mass[i] = mass[i-1] - dt*zs[i];
90 }
91
solve(V & c,V & u,const V & z,double & tol)92 void solve( V& c, V& u, const V& z, double& tol ) {
93
94 auto& us = getVector(u);
95
96 us[0] = -ve*std::log(mass[0]/mt) - gdt;
97 for( int i=1; i<N; ++i )
98 us[i] = us[i-1] - ve*std::log(mass[i]/mass[i-1]) - gdt;
99
100 value(c,u,z,tol);
101 }
102
value(V & c,const V & u,const V & z,double & tol)103 void value( V& c, const V& u, const V&z, double &tol ) {
104
105 auto& cs = getVector(c); auto& us = getVector(u);
106
107 cs[0] = us[0] + ve*std::log(mass[0]/mt) + gdt;
108
109 for( int i=1; i<N; ++i )
110 cs[i] = us[i] - us[i-1] + ve*std::log(mass[i]/mass[i-1]) + gdt;
111 }
112
applyJacobian_1(V & jv,const V & v,const V & u,const V & z,double & tol)113 void applyJacobian_1( V& jv, const V& v, const V& u, const V& z, double& tol ) {
114
115 auto& jvs = getVector(jv); auto& vs = getVector(v);
116
117 jvs[0] = vs[0];
118 for( int i=1; i<N; ++i ) jvs[i] = vs[i] - vs[i-1];
119 }
120
applyAdjointJacobian_1(V & ajv,const V & v,const V & u,const V & z,double & tol)121 void applyAdjointJacobian_1( V& ajv, const V& v, const V& u, const V& z, double& tol ) {
122
123 auto& ajvs = getVector(ajv); auto& vs = getVector(v);
124
125 ajvs[N-1] = vs[N-1];
126 for( int i=N-2; i>=0; --i ) ajvs[i] = vs[i] - vs[i+1];
127 }
128
applyInverseJacobian_1(V & ijv,const V & v,const V & u,const V & z,double & tol)129 void applyInverseJacobian_1( V& ijv, const V& v, const V& u, const V& z, double &tol ) {
130
131 auto& ijvs = getVector(ijv); auto& vs = getVector(v);
132
133 ijvs[0] = vs[0];
134 for( int i=1; i<N; ++i ) ijvs[i] = ijvs[i-1] + vs[i];
135 }
136
applyInverseAdjointJacobian_1(V & ijv,const V & v,const V & u,const V & z,double & tol)137 void applyInverseAdjointJacobian_1( V& ijv, const V& v, const V& u, const V& z, double &tol ) {
138
139 auto& ijvs = getVector(ijv); auto& vs = getVector(v);
140
141 ijvs[N-1] = vs[N-1];
142 for( int i=N-2; i>=0; --i ) ijvs[i] = ijvs[i+1] + vs[i];
143 }
144
applyJacobian_2(V & jv,const V & v,const V & u,const V & z,double & tol)145 void applyJacobian_2( V& jv, const V& v, const V& u, const V& z, double& tol ) {
146
147 auto& jvs = getVector(jv); auto& vs = getVector(v);
148
149 double q{-ve*dt*vs[0]};
150 jvs[0] = q/mass[0];
151 for( int i=1; i<N; ++i ) {
152 jvs[i] = -q/mass[i-1]; q -= ve*dt*vs[i];
153 jvs[i] += q/mass[i];
154 }
155 }
156
applyAdjointJacobian_2(V & ajv,const V & v,const V & u,const V & z,double & tol)157 void applyAdjointJacobian_2( V& ajv, const V& v, const V& u, const V& z, double& tol ) {
158 auto& ajvs = getVector(ajv); auto& vs = getVector(v);
159
160 ajvs[N-1] = -ve*dt*vs[N-1]/mass[N-1];
161
162 for( int i=N-2; i>=0; --i )
163 ajvs[i] = ajvs[i+1]-ve*dt*(vs[i]-vs[i+1])/mass[i];
164 }
165 };
166
167 } // namespace Rocket
168