PaGMO  1.1.5
sample_return.cpp
1 /*****************************************************************************
2  * Copyright (C) 2004-2015 The PaGMO development team, *
3  * Advanced Concepts Team (ACT), European Space Agency (ESA) *
4  * *
5  * https://github.com/esa/pagmo *
6  * *
7  * act@esa.int *
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  * This program is distributed in the hope that it will be useful, *
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
17  * GNU General Public License for more details. *
18  * *
19  * You should have received a copy of the GNU General Public License *
20  * along with this program; if not, write to the *
21  * Free Software Foundation, Inc., *
22  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
23  *****************************************************************************/
24 
25 #include <string>
26 
27 #include "sample_return.h"
28 #include <keplerian_toolbox/core_functions/convert_dates.h>
29 
30 
31 
32 namespace pagmo { namespace problem {
33 
34 static const int EA[2] = {3,10};
35 static const int AE[2] = {10,3};
36 
38 
41 sample_return::sample_return(const ::kep_toolbox::planet::base &asteroid, const double &Tmax):base(12), m_target(asteroid.clone()),
42  m_leg1(total_DV_rndv,EA,2,0,0,0,0,0), m_leg2(total_DV_rndv,AE,2,0,0,0,0,0),x_leg1(6),x_leg2(6),m_Tmax(Tmax)
43 {
44  ::kep_toolbox::epoch start_lw(2020,1,1);
45  ::kep_toolbox::epoch end_lw(2035,1,1);
46  const double lb[12] = {start_lw.mjd2000(),0 ,0,0,0.1 ,0.001, 20 ,0,0,0,0.1 ,0.001};
47  const double ub[12] = {end_lw.mjd2000() ,5.5,1,1,0.7, 0.999, 50 ,6,1,1,1 ,0.999};
48  set_bounds(lb,lb+12,ub,ub+12);
49 
50  for (int i = 0;i<6;++i)
51  {
52  m_leg1.asteroid.keplerian[i] = m_target->compute_elements(::kep_toolbox::epoch(0))[i];
53  m_leg2.asteroid.keplerian[i] = m_target->compute_elements(::kep_toolbox::epoch(0))[i];
54  }
55 
56  //convert to JPL unit format .... (check mgadsm)
57  m_leg1.asteroid.keplerian[0] /= ASTRO_AU;
58  m_leg2.asteroid.keplerian[0] /= ASTRO_AU;
59  for (int i = 2;i<6;++i)
60  {
61  m_leg1.asteroid.keplerian[i] *= ASTRO_RAD2DEG;
62  m_leg2.asteroid.keplerian[i] *= ASTRO_RAD2DEG;
63  }
64 
65  m_leg1.asteroid.epoch = ::kep_toolbox::mjd20002mjd(0);
66  m_leg2.asteroid.epoch = ::kep_toolbox::mjd20002mjd(0);
67 }
68 
71 {
72  return base_ptr(new sample_return(*this));
73 }
74 
77 {
78  //We split the decision vector in the two legs
79  std::copy(x.begin(),x.begin()+6,x_leg1.begin());
80  std::copy(x.begin()+6,x.begin()+12,x_leg2.begin());
81 
82  x_leg1[4] = x_leg1[4] * m_Tmax;
83  x_leg2[4] = (m_Tmax - x_leg1[4] - x_leg2[0]) * x_leg2[4];
84 
85  //We account for the waiting time
86  x_leg2[0] += x_leg1[0] + x_leg1[4];
87  double dummy = 0;
88  MGA_DSM(x_leg1, m_leg1, dummy);
89  MGA_DSM(x_leg2, m_leg2, dummy);
90  f[0] = m_leg1.DV[0] + m_leg1.DV[1] + m_leg1.DV[2] +
91  m_leg2.DV[0] + m_leg2.DV[1] + std::max(0.0,m_leg2.DV[2] - 5.5);
92 
93 }
94 
96 
105 std::string sample_return::pretty(const std::vector<double> &x) const
106 {
107  std::ostringstream s;
108  s.precision(15);
109  s << std::scientific;
110  //We split the decision vector in the two legs
111  std::copy(x.begin(),x.begin()+6,x_leg1.begin());
112  std::copy(x.begin()+6,x.begin()+12,x_leg2.begin());
113 
114  x_leg1[4] = x_leg1[4] * m_Tmax;
115  x_leg2[4] = (m_Tmax - x_leg1[4] - x_leg2[0]) * x_leg2[4];
116 
117  //We account for the waiting time
118  x_leg2[0] += x_leg1[0] + x_leg1[4];
119  double dummy = 0;
120  MGA_DSM(x_leg1, m_leg1, dummy);
121  MGA_DSM(x_leg2, m_leg2, dummy);
122 
123  s << "Departure epoch (mjd2000):\t" << x[0] << '\n';
124  s << "Departure epoch:\t\t" << ::kep_toolbox::epoch(x[0],::kep_toolbox::epoch::MJD2000) << '\n';
125  s << "Escape velocity:\t\t" << m_leg1.DV[0] << '\n';
126  s << "dsm1 epoch:\t\t\t" << ::kep_toolbox::epoch(x[0] + x[5]*x[4],::kep_toolbox::epoch::MJD2000) << '\n';
127  s << "dsm1 magnitude\t\t\t" << m_leg1.DV[1] << " \n";
128  s << "Asteroid arrival epoch: \t" << ::kep_toolbox::epoch(x[0] + x[4],::kep_toolbox::epoch::MJD2000) << '\n';
129  s << "Breaking manouvre:\t\t" << m_leg1.DV[2] << '\n' << std::endl;
130 
131  s << "Departure epoch (mjd2000):\t" << x[0] + x[4] + x[6] << '\n';
132  s << "Departure epoch:\t\t" << ::kep_toolbox::epoch(x[0] + x[4] + x[6],::kep_toolbox::epoch::MJD2000) << '\n';
133  s << "Escape velocity:\t\t" << m_leg2.DV[0] << '\n';
134  s << "dsm2 epoch:\t\t\t" << ::kep_toolbox::epoch(x[0]+x[4]+x[6]+x[11]*x[10],::kep_toolbox::epoch::MJD2000) << '\n';
135  s << "dsm2 magnitude\t\t\t" << m_leg2.DV[1] << " \n";
136  s << "Earth arrival epoch: \t\t" << ::kep_toolbox::epoch(x[0]+x[4]+x[6]+x[10],::kep_toolbox::epoch::MJD2000) << '\n';
137  s << "Arrival Vinf:\t\t\t" << m_leg2.DV[2] << '\n';
138  s << "Total time of flight:\t\t" << x[4]+x[6]+x[10] << '\n' << std::endl;
139 
140  s << "Earth-ephemerides at departure:\t\t" << m_leg1.r[0][0] << " " << m_leg1.r[0][1] << " " << m_leg1.r[0][2] << std::endl;
141  s << "Asteroid-ephemerides at arrival:\t" << m_leg1.r[1][0] << " " << m_leg1.r[1][1] << " " << m_leg1.r[1][2] << std::endl;
142  s << "Asteroid-ephemerides at departure:\t" << m_leg2.r[0][0] << " " << m_leg2.r[0][1] << " " << m_leg2.r[0][2] << std::endl;
143  s << "Earth-ephemerides at arrival:\t\t" << m_leg2.r[1][0] << " " << m_leg2.r[1][1] << " " << m_leg2.r[1][2] << std::endl;
144 
145 
146 
147  return s.str();
148 }
149 
151 std::vector<double> sample_return::get_delta_v(const std::vector<double> &x) const {
152  //We split the decision vector in the two legs
153  std::copy(x.begin(),x.begin()+6,x_leg1.begin());
154  std::copy(x.begin()+6,x.begin()+12,x_leg2.begin());
155 
156  x_leg1[4] = x_leg1[4] * m_Tmax;
157  x_leg2[4] = (m_Tmax - x_leg1[4] - x_leg2[0]) * x_leg2[4];
158 
159  //We account for the waiting time
160  x_leg2[0] += x_leg1[0] + x_leg1[4];
161  double dummy = 0;
162  MGA_DSM(x_leg1, m_leg1, dummy);
163  MGA_DSM(x_leg2, m_leg2, dummy);
164  std::vector<double> retval;
165  for (int i=0;i<3;++i) retval.push_back(m_leg1.DV[i]);
166  for (int i=0;i<3;++i) retval.push_back(m_leg2.DV[i]);
167  return retval;
168 }
169 
171 
174 void sample_return::set_sparsity(int &lenG, std::vector<int> &iGfun, std::vector<int> &jGvar) const
175 {
176  lenG=12;
177  iGfun.resize(lenG);
178  jGvar.resize(lenG);
179  for (int i = 0; i<lenG; ++i)
180  {
181  iGfun[i] = 0;
182  jGvar[i] = i;
183  }
184 }
185 
186 std::string sample_return::get_name() const
187 {
188  return "Sample return";
189 }
190 
191 }} //namespaces
192 
193 BOOST_CLASS_EXPORT_IMPLEMENT(pagmo::problem::sample_return)
Root PaGMO namespace.
boost::shared_ptr< base > base_ptr
Alias for shared pointer to base problem.
Definition: problem/base.h:62
std::vector< double > decision_vector
Decision vector type.
Definition: types.h:40
Base problem class.
Definition: problem/base.h:148
Human mission to asteroids.
Definition: sample_return.h:50
std::vector< double > get_delta_v(const std::vector< double > &x) const
Computes the Delta-Vs.
sample_return(const ::kep_toolbox::planet::base &asteroid=::kep_toolbox::planet::mpcorb(), const double &Tmax=600)
Constructor.
std::vector< double > fitness_vector
Fitness vector type.
Definition: types.h:42
void objfun_impl(fitness_vector &, const decision_vector &) const
Implementation of the objective function.
std::string get_name() const
Get problem's name.
void set_sparsity(int &, std::vector< int > &, std::vector< int > &) const
Implementation of the sparsity structure.
base_ptr clone() const
Clone method.
void set_bounds(const decision_vector &, const decision_vector &)
Bounds setter from pagmo::decision_vector.
std::string pretty(const std::vector< double > &x) const
Outputs a stream with the trajectory data.