29 #include <boost/math/constants/constants.hpp>
30 #include <boost/array.hpp>
31 #include <keplerian_toolbox/core_functions/propagate_lagrangian.h>
32 #include <keplerian_toolbox/core_functions/fb_prop.h>
33 #include <keplerian_toolbox/lambert_problem.h>
35 #include "mga_1dsm_tof.h"
38 namespace pagmo {
namespace problem {
59 const kep_toolbox::epoch t0_l,
const kep_toolbox::epoch t0_u,
60 const std::vector<boost::array<double,2> > tof,
61 const double vinf_l,
const double vinf_u,
62 const bool mo,
const bool add_vinf_dep,
const bool add_vinf_arr) :
63 base(6 + (int)(seq.size()-2) * 4, 0, 1 + (int)mo,0,0,0.0), m_n_legs(seq.size()-1), m_add_vinf_dep(add_vinf_dep), m_add_vinf_arr(add_vinf_arr)
66 std::vector<double> mus(seq.size());
67 for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i< seq.size(); ++i) {
68 mus[i] = seq[i]->get_mu_central_body();
70 if ((
unsigned int)std::count(mus.begin(), mus.end(), mus[0]) != mus.size()) {
71 pagmo_throw(value_error,
"The planets do not all have the same mu_central_body");
74 if (seq.size()!= (tof.size()+1)) {
75 pagmo_throw(value_error,
"The number of time of flights is not compatible with the sequence length");
79 for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i < seq.size(); ++i) {
80 m_seq.push_back(seq[i]->
clone());
86 lb[0] = t0_l.mjd2000(); ub[0] = t0_u.mjd2000();
89 lb[1] = 0; lb[2] = 0; ub[1] = 1; ub[2] = 1;
90 lb[3] = vinf_l * 1000; ub[3] = vinf_u * 1000;
91 lb[4] = 1e-5; ub[4] = 1-1e-5;
92 lb[5] = tof[0][0]; ub[5] = tof[0][1];
95 for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i < m_n_legs - 1; ++i) {
96 lb[6+4*i] = - 2 * boost::math::constants::pi<double>(); ub[6+4*i] = 2 * boost::math::constants::pi<double>();
97 lb[7+4*i] = 1.1; ub[7+4*i] = 100;
98 lb[8+4*i] = 1e-5; ub[8+4*i] = 1-1e-5;
99 lb[9+4*i] = tof[i+1][0]; ub[9+4*i] = tof[i+1][1];
102 for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 1; i < m_n_legs; ++i) {
103 lb[3 + 4*i] = m_seq[i]->get_safe_radius() / m_seq[i]->get_radius();
113 for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i < p.m_seq.size();++i) {
114 m_seq.push_back(p.m_seq[i]->clone());
129 double common_mu = m_seq[0]->get_mu_central_body();
131 std::vector<double> T(m_n_legs,0.0);
133 for (
size_t i = 0; i < m_n_legs; ++i) {
138 std::vector<kep_toolbox::epoch> t_P(m_n_legs + 1);
139 std::vector<kep_toolbox::array3D> r_P(m_n_legs + 1);
140 std::vector<kep_toolbox::array3D> v_P(m_n_legs + 1);
141 std::vector<double> DV(m_n_legs + 1);
142 for (
size_t i = 0; i<(m_n_legs + 1); ++i) {
143 t_P[i] = kep_toolbox::epoch(x[0] + std::accumulate(T.begin(), T.begin()+i, 0.0));
144 m_seq[i]->eph(t_P[i], r_P[i], v_P[i]);
148 double theta = 2*boost::math::constants::pi<double>()*x[1];
149 double phi = acos(2*x[2]-1)-boost::math::constants::pi<double>() / 2;
151 kep_toolbox::array3D Vinf = { {x[3]*cos(phi)*cos(theta), x[3]*cos(phi)*sin(theta), x[3]*sin(phi)} };
152 kep_toolbox::array3D v0;
153 kep_toolbox::sum(v0, v_P[0], Vinf);
154 kep_toolbox::array3D r(r_P[0]), v(v0);
155 kep_toolbox::propagate_lagrangian(r,v,x[4]*T[0]*ASTRO_DAY2SEC,common_mu);
158 double dt = (1-x[4])*T[0]*ASTRO_DAY2SEC;
159 kep_toolbox::lambert_problem l(r,r_P[1],dt,common_mu);
160 kep_toolbox::array3D v_end_l = l.get_v2()[0];
161 kep_toolbox::array3D v_beg_l = l.get_v1()[0];
164 kep_toolbox::diff(v, v_beg_l, v);
165 DV[0] = kep_toolbox::norm(v);
168 kep_toolbox::array3D v_out;
169 for (
size_t i = 1; i<m_n_legs; ++i) {
171 kep_toolbox::fb_prop(v_out, v_end_l, v_P[i], x[7+(i-1)*4] * m_seq[i]->get_radius(), x[6+(i-1)*4], m_seq[i]->get_mu_self());
175 kep_toolbox::propagate_lagrangian(r,v,x[8+(i-1)*4]*T[i]*ASTRO_DAY2SEC,common_mu);
178 dt = (1-x[8+(i-1)*4])*T[i]*ASTRO_DAY2SEC;
179 kep_toolbox::lambert_problem l2(r,r_P[i+1],dt,common_mu);
180 v_end_l = l2.get_v2()[0];
181 v_beg_l = l2.get_v1()[0];
184 kep_toolbox::diff(v, v_beg_l, v);
185 DV[i] = kep_toolbox::norm(v);
189 kep_toolbox::diff(v, v_end_l, v_P[m_n_legs]);
190 DV[m_n_legs] = kep_toolbox::norm(v);
194 f[0] = std::accumulate(DV.begin(),DV.end()-1,0.0);
195 if (m_add_vinf_dep) {
198 if (m_add_vinf_arr) {
199 f[0] += DV[DV.size()-1];
202 f[1] = std::accumulate(T.begin(),T.end(),0.0);
206 f[0] = boost::numeric::bounds<double>::highest();
208 f[1] = boost::numeric::bounds<double>::highest();
226 double common_mu = m_seq[0]->get_mu_central_body();
228 std::ostringstream s;
230 s << std::scientific;
233 std::vector<double> T(m_n_legs,0.0);
235 for (
size_t i = 0; i < m_n_legs; ++i) {
240 std::vector<kep_toolbox::epoch> t_P(m_n_legs + 1);
241 std::vector<kep_toolbox::array3D> r_P(m_n_legs + 1);
242 std::vector<kep_toolbox::array3D> v_P(m_n_legs + 1);
243 std::vector<double> DV(m_n_legs + 1);
244 for (
size_t i = 0; i<(m_n_legs + 1); ++i) {
245 t_P[i] = kep_toolbox::epoch(x[0] + std::accumulate(T.begin(), T.begin()+i, 0.0));
246 m_seq[i]->eph(t_P[i], r_P[i], v_P[i]);
250 s <<
"First Leg: " << m_seq[0]->get_name() <<
" to " << m_seq[1]->get_name() << std::endl;
251 s <<
"Departure: " << t_P[0] <<
" (" << t_P[0].mjd2000() <<
" mjd2000)" << std::endl;
252 s <<
"Duration: " << T[0] <<
" days" << std::endl;
253 double theta = 2*boost::math::constants::pi<double>()*x[1];
254 double phi = acos(2*x[2]-1)-boost::math::constants::pi<double>() / 2;
256 kep_toolbox::array3D Vinf = { {x[3]*cos(phi)*cos(theta), x[3]*cos(phi)*sin(theta), x[3]*sin(phi)} };
258 kep_toolbox::array3D v_dep, v_misc;
259 kep_toolbox::sum(v_dep, v_P[0], Vinf);
260 kep_toolbox::array3D r(r_P[0]), v(v_dep);
262 if(extended_output) s <<
"r_dep: " << r <<
"\nv_dep: " << v << std::endl;
263 s <<
"VINF_dep: " << x[3] <<
" m/sec" << std::endl;
264 kep_toolbox::propagate_lagrangian(r,v,x[4]*T[0]*ASTRO_DAY2SEC,common_mu);
267 double dt = (1-x[4])*T[0]*ASTRO_DAY2SEC;
268 kep_toolbox::lambert_problem l(r,r_P[1],dt,common_mu);
269 kep_toolbox::array3D v_end_l = l.get_v2()[0];
270 kep_toolbox::array3D v_beg_l = l.get_v1()[0];
272 if(extended_output) s <<
"r_arr: " << r_P[1] <<
"\nv_arr: " << v_end_l << std::endl;
273 kep_toolbox::diff(v_misc, v_end_l,v_P[1]);
274 if(extended_output) s <<
"VINF_arr: " << kep_toolbox::norm(v_misc) <<
" m/sec" << std::endl;
277 s <<
"DSM after " << x[4]*T[0] <<
" days" << std::endl;
278 kep_toolbox::diff(v, v_beg_l, v);
279 DV[0] = kep_toolbox::norm(v);
280 s <<
"DSM magnitude: " << DV[0] <<
" m/s" << std::endl;
283 kep_toolbox::array3D v_out;
284 for (
size_t i = 1; i<m_n_legs; ++i) {
285 s <<
"\nleg no: " << i+1 <<
": " << m_seq[i]->get_name() <<
" to " << m_seq[i+1]->get_name() << std::endl;
286 s <<
"Duration: " << T[i] <<
" days" << std::endl;
287 s <<
"Fly-by epoch: " << t_P[i] <<
" (" << t_P[i].mjd2000() <<
" mjd2000) " << std::endl;
288 s <<
"Fly-by radius: " << x[7+(i-1)*4] <<
" planetary radii" << std::endl;
291 kep_toolbox::fb_prop(v_out, v_end_l, v_P[i], x[7+(i-1)*4] * m_seq[i]->get_radius(), x[6+(i-1)*4], m_seq[i]->get_mu_self());
295 if(extended_output) {
296 s <<
"r_dep: " << r_P[i] <<
"\nv_dep: " << v_out << std::endl;
297 kep_toolbox::diff(v_misc, v_out,v_P[i]);
298 s <<
"VINF_dep: " << kep_toolbox::norm(v_misc) <<
" m/sec" << std::endl;
300 kep_toolbox::propagate_lagrangian(r,v,x[8+(i-1)*4]*T[i]*ASTRO_DAY2SEC,common_mu);
303 dt = (1-x[8+(i-1)*4])*T[i]*ASTRO_DAY2SEC;
304 kep_toolbox::lambert_problem l2(r,r_P[i+1],dt,common_mu);
305 v_end_l = l2.get_v2()[0];
306 v_beg_l = l2.get_v1()[0];
309 kep_toolbox::diff(v, v_beg_l, v);
310 DV[i] = kep_toolbox::norm(v);
311 if(extended_output) {
312 s <<
"r_arr: " << r_P[i+1] <<
"\nv_arr: " << v_end_l << std::endl;
313 kep_toolbox::diff(v_misc, v_end_l,v_P[i+1]);
314 s <<
"VINF_arr: " << kep_toolbox::norm(v_misc) <<
" m/sec" << std::endl;
316 s <<
"DSM after " << x[8+(i-1)*4]*T[i] <<
" days" << std::endl;
317 s <<
"DSM magnitude: " << DV[i] <<
"m/s" << std::endl;
321 kep_toolbox::diff(v, v_end_l, v_P[m_n_legs]);
322 DV[m_n_legs] = kep_toolbox::norm(v);
323 s <<
"\nArrival at " << m_seq[m_n_legs]->get_name() << std::endl;
324 s <<
"Arrival epoch: " << t_P[m_n_legs] <<
" (" << t_P[m_n_legs].mjd2000() <<
" mjd2000) " << std::endl;
325 s <<
"Arrival Vinf: " << DV[m_n_legs] <<
"m/s" << std::endl;
326 s <<
"Total mission time: " << std::accumulate(T.begin(),T.end(),0.0)/365.25 <<
" years" << std::endl;
331 return "MGA-1DSM (tof-encoding)";
341 if (m_seq.size()!= (tof.size()+1)) {
342 pagmo_throw(value_error,
"The size of the time of flight is inconsistent");
345 for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i < m_n_legs; ++i) {
376 std::vector<std::vector<double> > retval;
377 std::vector<double> tmp;
380 for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i < m_n_legs; ++i) {
381 tmp.push_back(lb[5+(i*4)]);
382 tmp.push_back(ub[5+(i*4)]);
383 retval.push_back(tmp);
404 std::ostringstream oss;
405 oss <<
"\n\tSequence: ";
406 for (
size_t i = 0; i<m_seq.size() ;++i) {
407 oss << m_seq[i]->get_name() <<
" ";
409 oss <<
"\n\tAdd launcher vinf to the objective?: " << (m_add_vinf_dep?
" True":
" False") << std::endl;
410 oss <<
"\n\tAdd arrival vinf to the objective?: " << (m_add_vinf_arr?
" True":
" False") << std::endl;
boost::shared_ptr< base > base_ptr
Alias for shared pointer to base problem.
std::vector< double > decision_vector
Decision vector type.
mga_1dsm_tof(const std::vector< kep_toolbox::planet::planet_ptr >=construct_default_sequence(), const kep_toolbox::epoch t0_l=kep_toolbox::epoch(0), const kep_toolbox::epoch t0_r=kep_toolbox::epoch(1000), const std::vector< boost::array< double, 2 > >=construct_default_tof(), const double vinf_l=0.5, const double vinf_u=2.5, const bool mo=false, const bool add_vinf_dep=false, const bool add_vinf_arr=true)
Constructor.
std::vector< kep_toolbox::planet::planet_ptr > get_sequence() const
Gets the planetary sequence defining the interplanetary mga-1dsm mission.
std::string human_readable_extra() const
Extra human readable info for the problem.
std::string get_name() const
Get problem's name.
base_ptr clone() const
Clone method.
void set_launch_window(const kep_toolbox::epoch &, const kep_toolbox::epoch &)
Sets the mission launch window.
void set_tof(const std::vector< boost::array< double, 2 > >)
Sets the mission time of flight.
A generic MGA-1DSM Problem (tof encoding)
std::string pretty(const std::vector< double > &x, bool extended_output=false) const
Outputs a stream with the trajectory data.
void set_ub(const decision_vector &)
Set upper bounds from pagmo::decision_vector.
std::vector< double > fitness_vector
Fitness vector type.
void objfun_impl(fitness_vector &, const decision_vector &) const
Implementation of the objective function.
void set_vinf(const double)
Sets the launch hyperbolic velocity.
const decision_vector & get_ub() const
Upper bounds getter.
std::vector< std::vector< double > > get_tof() const
Gets the sequence of time of flights for the mga-1dsm mission.
f_size_type get_f_dimension() const
Return fitness dimension.
const decision_vector & get_lb() const
Lower bounds getter.
void set_bounds(const decision_vector &, const decision_vector &)
Bounds setter from pagmo::decision_vector.
decision_vector::size_type size_type
Problem's size type: the same as pagmo::decision_vector's size type.