PaGMO  1.1.5
mga_1dsm_tof.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 #include <vector>
27 #include <numeric>
28 #include <cmath>
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>
34 
35 #include "mga_1dsm_tof.h"
36 
37 
38 namespace pagmo { namespace problem {
39 
40 
42 
58 mga_1dsm_tof::mga_1dsm_tof(const std::vector<kep_toolbox::planet::planet_ptr> seq,
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)
64 {
65  // We check that all planets have equal central body
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();
69  }
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");
72  }
73  // And that the time of flights have the correct dimension
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");
76  }
77 
78  // Filling in the planetary sequence data member. This requires to construct the polymorphic planets via their clone method
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());
81  }
82  // Now setting the problem bounds
83  size_type dim(6 + (m_n_legs-1) * 4);
84  decision_vector lb(dim), ub(dim);
85  // First leg
86  lb[0] = t0_l.mjd2000(); ub[0] = t0_u.mjd2000();
87 
88  // TODO: implement
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];
93 
94  // Successive legs
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];
100  }
101  // Adjusting the minimum allowed fly-by rp to the one defined in the kep_toolbox::planet class
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();
104  }
105 
106  set_bounds(lb,ub);
107 
108 }
109 
111 mga_1dsm_tof::mga_1dsm_tof(const mga_1dsm_tof &p) : base(p.get_dimension(), 0, p.get_f_dimension(),0,0,0.0), m_seq(), m_n_legs(p.m_n_legs), m_add_vinf_dep(p.m_add_vinf_dep), m_add_vinf_arr(p.m_add_vinf_arr)
112 {
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());
115  }
116  set_bounds(p.get_lb(),p.get_ub());
117 }
118 
121 {
122  return base_ptr(new mga_1dsm_tof(*this));
123 }
124 
127 {
128 try {
129  double common_mu = m_seq[0]->get_mu_central_body();
130  // 1 - we 'decode' the chromosome recording the various times of flight (days) in the list T
131  std::vector<double> T(m_n_legs,0.0);
132 
133  for (size_t i = 0; i < m_n_legs; ++i) {
134  T[i] = x[5 + i*4];
135  }
136 
137  // 2 - We compute the epochs and ephemerides of the planetary encounters
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]);
145  }
146 
147  // 3 - We start with the first leg
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;
150 
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);
156 
157  // Lambert arc to reach seq[1]
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];
162 
163  // First DSM occuring at time nu1*T1
164  kep_toolbox::diff(v, v_beg_l, v);
165  DV[0] = kep_toolbox::norm(v);
166 
167  // 4 - And we proceed with each successive leg (if any)
168  kep_toolbox::array3D v_out;
169  for (size_t i = 1; i<m_n_legs; ++i) {
170  // Fly-by
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());
172  // s/c propagation before the DSM
173  r = r_P[i];
174  v = v_out;
175  kep_toolbox::propagate_lagrangian(r,v,x[8+(i-1)*4]*T[i]*ASTRO_DAY2SEC,common_mu);
176 
177  // Lambert arc to reach Earth during (1-nu2)*T2 (second segment)
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];
182 
183  // DSM occuring at time nu2*T2
184  kep_toolbox::diff(v, v_beg_l, v);
185  DV[i] = kep_toolbox::norm(v);
186  }
187 
188  // Last Delta-v
189  kep_toolbox::diff(v, v_end_l, v_P[m_n_legs]);
190  DV[m_n_legs] = kep_toolbox::norm(v);
191 
192 
193  // Now we return the objective(s) function
194  f[0] = std::accumulate(DV.begin(),DV.end()-1,0.0);
195  if (m_add_vinf_dep) {
196  f[0] += x[3];
197  }
198  if (m_add_vinf_arr) {
199  f[0] += DV[DV.size()-1];
200  }
201  if (get_f_dimension() == 2){
202  f[1] = std::accumulate(T.begin(),T.end(),0.0);
203  }
204 //Here the lambert solver or the lagrangian propagator went wrong
205 } catch (...) {
206  f[0] = boost::numeric::bounds<double>::highest();
207  if (get_f_dimension() == 2){
208  f[1] = boost::numeric::bounds<double>::highest();
209  }
210 }
211 }
212 
214 
224 std::string mga_1dsm_tof::pretty(const std::vector<double> &x, bool extended_output) const {
225 
226  double common_mu = m_seq[0]->get_mu_central_body();
227  // We set the std output format
228  std::ostringstream s;
229  s.precision(15);
230  s << std::scientific;
231 
232  // 1 - we 'decode' the chromosome recording the various times of flight (days) in the list T
233  std::vector<double> T(m_n_legs,0.0);
234 
235  for (size_t i = 0; i < m_n_legs; ++i) {
236  T[i] = x[5 + i*4];
237  }
238 
239  // 2 - We compute the epochs and ephemerides of the planetary encounters
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]);
247  }
248 
249  // 3 - We start with the first leg
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;
255 
256  kep_toolbox::array3D Vinf = { {x[3]*cos(phi)*cos(theta), x[3]*cos(phi)*sin(theta), x[3]*sin(phi)} };
257 
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);
261 
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);
265 
266  // Lambert arc to reach seq[1]
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];
271 
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;
275 
276  // First DSM occuring at time nu1*T1
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;
281 
282  // 4 - And we proceed with each successive leg (if any)
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;
289 
290  // Fly-by
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());
292  // s/c propagation before the DSM
293  r = r_P[i];
294  v = v_out;
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;
299  }
300  kep_toolbox::propagate_lagrangian(r,v,x[8+(i-1)*4]*T[i]*ASTRO_DAY2SEC,common_mu);
301 
302  // Lambert arc to reach Earth during (1-nu2)*T2 (second segment)
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];
307 
308  // DSM occuring at time nu2*T2
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;
315  }
316  s << "DSM after " << x[8+(i-1)*4]*T[i] << " days" << std::endl;
317  s << "DSM magnitude: " << DV[i] << "m/s" << std::endl;
318  }
319 
320  // Last Delta-v
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;
327  return s.str();
328 }
329 std::string mga_1dsm_tof::get_name() const
330 {
331  return "MGA-1DSM (tof-encoding)";
332 }
333 
335 
340 void mga_1dsm_tof::set_tof(const std::vector<boost::array<double,2> > tof) {
341  if (m_seq.size()!= (tof.size()+1)) {
342  pagmo_throw(value_error,"The size of the time of flight is inconsistent");
343  }
344  // setting bounds
345  for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i < m_n_legs; ++i) {
346  set_bounds(5 + i*4,tof[i][0],tof[i][1]);
347  }
348 }
349 
351 
357 void mga_1dsm_tof::set_launch_window(const kep_toolbox::epoch& start, const kep_toolbox::epoch& end) {
358  set_bounds(0,start.mjd2000(),end.mjd2000());
359 }
360 
362 
367 void mga_1dsm_tof::set_vinf(const double vinf) {
368  set_ub(3,vinf*1000);
369 }
370 
372 
375 std::vector<std::vector<double> > mga_1dsm_tof::get_tof() const {
376  std::vector<std::vector<double> > retval;
377  std::vector<double> tmp;
378  const decision_vector lb = get_lb();
379  const decision_vector ub = get_ub();
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);
384  tmp.clear();
385  }
386  return retval;
387 }
388 
390 
393 std::vector<kep_toolbox::planet::planet_ptr> mga_1dsm_tof::get_sequence() const {
394  return m_seq;
395 }
396 
398 
403 {
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() << " ";
408  }
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;
411  return oss.str();
412 }
413 
414 }} //namespaces
415 
416 BOOST_CLASS_EXPORT_IMPLEMENT(pagmo::problem::mga_1dsm_tof)
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
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.
Base problem class.
Definition: problem/base.h:148
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)
Definition: mga_1dsm_tof.h:63
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.
Definition: types.h:42
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.
Definition: problem/base.h:160