PaGMO  1.1.5
mga_incipit_cstrs.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 <boost/math/constants/constants.hpp>
27 #include <vector>
28 #include <numeric>
29 #include <cmath>
30 #include <keplerian_toolbox/core_functions/propagate_lagrangian.h>
31 #include <keplerian_toolbox/core_functions/fb_prop.h>
32 #include <keplerian_toolbox/core_functions/closest_distance.h>
33 #include <keplerian_toolbox/lambert_problem.h>
34 
35 #include "mga_incipit_cstrs.h"
36 
37 
38 namespace pagmo { namespace problem {
39 
40 
42 
57  const std::vector<kep_toolbox::planet::planet_ptr> seq,
58  const kep_toolbox::epoch t0_l,
59  const kep_toolbox::epoch t0_u,
60  const std::vector<std::vector<double> > tof,
61  double Tmax,
62  double Dmin) : base(4*seq.size(),0,1,2,2,1E-3), m_tof(tof), m_tmax(Tmax), m_dmin(Dmin)
63 {
64  // We check that all planets have equal central body
65  std::vector<double> mus(seq.size());
66  for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i< seq.size(); ++i) {
67  mus[i] = seq[i]->get_mu_central_body();
68  }
69  if ((unsigned int)std::count(mus.begin(), mus.end(), mus[0]) != mus.size()) {
70  pagmo_throw(value_error,"The planets do not all have the same mu_central_body");
71  }
72  // We check the consistency of the time of flights
73  if (tof.size() != seq.size()) {
74  pagmo_throw(value_error,"The time-of-flight vector (tof) has the wrong length");
75  }
76  for (size_t i = 0; i < tof.size(); ++i) {
77  if (tof[i].size()!=2) pagmo_throw(value_error,"Each element of the time-of-flight vector (tof) needs to have dimension 2 (lower and upper bound)");
78  }
79 
80  // Filling in the planetary sequence data member. This requires to construct the polymorphic planets via their clone method
81  for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i < seq.size(); ++i) {
82  m_seq.push_back(seq[i]->clone());
83  }
84 
85  // Now setting the problem bounds
86  size_type dim(4*m_tof.size());
87  decision_vector lb(dim), ub(dim);
88 
89  // First leg
90  lb[0] = t0_l.mjd2000(); ub[0] = t0_u.mjd2000();
91  lb[1] = 0; lb[2] = 0; ub[1] = 1; ub[2] = 1;
92  lb[3] = m_tof[0][0]; ub[3] = m_tof[0][1];
93 
94  // Successive legs
95  for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 1; i < m_tof.size(); ++i) {
96  lb[4*i] = - 2 * boost::math::constants::pi<double>(); ub[4*i] = 2 * boost::math::constants::pi<double>();
97  lb[4*i+1] = 1.1; ub[4*i+1] = 30;
98  lb[4*i+2] = 1e-5; ub[4*i+2] = 1-1e-5;
99  lb[4*i+3] = m_tof[i][0]; ub[4*i+3] = m_tof[i][1];
100  }
101 
102  // Adjusting the minimum and maximum allowed fly-by rp to the one defined in the kep_toolbox::planet class
103  for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i < m_tof.size()-1; ++i) {
104  lb[4*i+5] = m_seq[i]->get_safe_radius() / m_seq[i]->get_radius();
105  ub[4*i+5] = (m_seq[i]->get_radius() + 2000000) / m_seq[i]->get_radius(); //from gtoc6 problem description
106  }
107  set_bounds(lb,ub);
108 }
109 
112  base(p.get_dimension(),p.get_i_dimension(),p.get_f_dimension(),p.get_c_dimension(),p.get_ic_dimension(),p.get_c_tol()),
113  m_tof(p.m_tof),
114  m_tmax(p.m_tmax),
115  m_dmin(p.m_dmin)
116 {
117  for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i < p.m_seq.size();++i) {
118  m_seq.push_back(p.m_seq[i]->clone());
119  }
120  set_bounds(p.get_lb(),p.get_ub());
121 }
122 
125 {
126  return base_ptr(new mga_incipit_cstrs(*this));
127 }
128 
131 {
132 try {
133  double common_mu = m_seq[0]->get_mu_central_body();
134  // 1 - we 'decode' the chromosome recording the various times of flight (days) in the list T
135  std::vector<double> T(m_seq.size(),0.0);
136 
137  for (size_t i = 0; i<m_seq.size(); ++i) {
138  T[i] = x[4*i+3];
139  }
140  // 2 - We compute the epochs and ephemerides of the planetary encounters
141  std::vector<kep_toolbox::epoch> t_P(m_seq.size());
142  std::vector<kep_toolbox::array3D> r_P(m_seq.size());
143  std::vector<kep_toolbox::array3D> v_P(m_seq.size());
144  std::vector<double> DV(m_seq.size());
145  for (size_t i = 0; i<r_P.size(); ++i) {
146  t_P[i] = kep_toolbox::epoch(x[0] + std::accumulate(T.begin(),T.begin()+1+i,0.0));
147  m_seq[i]->eph(t_P[i], r_P[i], v_P[i]);
148  }
149 
150  // 3 - We start with the first leg
151  double theta = 2*boost::math::constants::pi<double>()*x[1];
152  double phi = acos(2*x[2]-1)-boost::math::constants::pi<double>() / 2;
153  double d,d2,ra,ra2;
154  kep_toolbox::array3D r = { {ASTRO_JR*1000*cos(phi)*sin(theta), ASTRO_JR*1000*cos(phi)*cos(theta), ASTRO_JR*1000*sin(phi)} };
155  kep_toolbox::array3D v;
156  kep_toolbox::lambert_problem l(r,r_P[0],T[0]*ASTRO_DAY2SEC,common_mu,false,false);
157  kep_toolbox::array3D v_beg_l = l.get_v1()[0];
158  kep_toolbox::array3D v_end_l = l.get_v2()[0];
159 
160  DV[0] = std::abs(kep_toolbox::norm(v_beg_l)-3400.0);
161 
162  // 4 - And we proceed with each successive leg (if any)
163  kep_toolbox::array3D v_out;
164  for (size_t i = 1; i<m_seq.size(); ++i) {
165  // Fly-by
166  kep_toolbox::fb_prop(v_out, v_end_l, v_P[i-1], x[4*i+1] * m_seq[i-1]->get_radius(), x[4*i], m_seq[i-1]->get_mu_self());
167  r = r_P[i-1];
168  v = v_out;
169  // s/c propagation before the DSM
170  kep_toolbox::propagate_lagrangian(r,v,x[4*i+2]*T[i]*ASTRO_DAY2SEC,common_mu);
171  kep_toolbox::closest_distance(d, ra, r_P[i-1], v_out, r, v, common_mu);
172 
173  // Lambert arc to reach Earth during (1-nu2)*T2 (second segment)
174  double dt = (1-x[4*i+2])*T[i]*ASTRO_DAY2SEC;
175  kep_toolbox::lambert_problem l2(r,r_P[i],dt,common_mu,false,false);
176  v_end_l = l2.get_v2()[0];
177  v_beg_l = l2.get_v1()[0];
178  kep_toolbox::closest_distance(d2,ra2,r,v_beg_l, r_P[i], v_end_l, common_mu);
179  if (d < d2)
180  {
181  d = d/ASTRO_JR;
182  } else {
183  d = d2/ASTRO_JR;
184  }
185 
186  // DSM occuring at time nu2*T2
187  kep_toolbox::diff(v_out, v_beg_l, v);
188  DV[i] = kep_toolbox::norm(v_out) + std::max((2.0-d),0.0) * 1000.0;
189  }
190  // Now we return the objective(s) function
191  f[0] = std::accumulate(DV.begin(),DV.end(),0.0);
192 //Here the lambert solver or the lagrangian propagator went wrong
193 } catch (...) {
194  f[0] = boost::numeric::bounds<double>::highest();
195 }
196 }
197 
200 {
201 try {
202  double common_mu = m_seq[0]->get_mu_central_body();
203  // 1 - we 'decode' the chromosome recording the various times of flight (days) in the list T
204  std::vector<double> T(m_seq.size(),0.0);
205 
206  for (size_t i = 0; i<m_seq.size(); ++i) {
207  T[i] = x[4*i+3];
208  }
209  // 2 - We compute the epochs and ephemerides of the planetary encounters
210  std::vector<kep_toolbox::epoch> t_P(m_seq.size());
211  std::vector<kep_toolbox::array3D> r_P(m_seq.size());
212  std::vector<kep_toolbox::array3D> v_P(m_seq.size());
213  std::vector<double> DV(m_seq.size());
214  for (size_t i = 0; i<r_P.size(); ++i) {
215  t_P[i] = kep_toolbox::epoch(x[0] + std::accumulate(T.begin(),T.begin()+1+i,0.0));
216  m_seq[i]->eph(t_P[i], r_P[i], v_P[i]);
217  }
218 
219  // 3 - We start with the first leg
220  double theta = 2*boost::math::constants::pi<double>()*x[1];
221  double phi = acos(2*x[2]-1)-boost::math::constants::pi<double>() / 2;
222  double d,d2,ra,ra2;
223  kep_toolbox::array3D r = { {ASTRO_JR*1000*cos(phi)*sin(theta), ASTRO_JR*1000*cos(phi)*cos(theta), ASTRO_JR*1000*sin(phi)} };
224  kep_toolbox::array3D v;
225  kep_toolbox::lambert_problem l(r,r_P[0],T[0]*ASTRO_DAY2SEC,common_mu,false,false);
226  kep_toolbox::array3D v_beg_l = l.get_v1()[0];
227  kep_toolbox::array3D v_end_l = l.get_v2()[0];
228 
229  DV[0] = std::abs(kep_toolbox::norm(v_beg_l)-3400.0);
230 
231  // 4 - And we proceed with each successive leg (if any)
232  kep_toolbox::array3D v_out;
233  for (size_t i = 1; i<m_seq.size(); ++i) {
234  // Fly-by
235  kep_toolbox::fb_prop(v_out, v_end_l, v_P[i-1], x[4*i+1] * m_seq[i-1]->get_radius(), x[4*i], m_seq[i-1]->get_mu_self());
236  r = r_P[i-1];
237  v = v_out;
238  // s/c propagation before the DSM
239  kep_toolbox::propagate_lagrangian(r,v,x[4*i+2]*T[i]*ASTRO_DAY2SEC,common_mu);
240  kep_toolbox::closest_distance(d, ra, r_P[i-1], v_out, r, v, common_mu);
241 
242  // Lambert arc to reach Earth during (1-nu2)*T2 (second segment)
243  double dt = (1-x[4*i+2])*T[i]*ASTRO_DAY2SEC;
244  kep_toolbox::lambert_problem l2(r,r_P[i],dt,common_mu,false,false);
245  v_end_l = l2.get_v2()[0];
246  v_beg_l = l2.get_v1()[0];
247  kep_toolbox::closest_distance(d2,ra2,r,v_beg_l, r_P[i], v_end_l, common_mu);
248  if (d < d2)
249  {
250  d = d/ASTRO_JR;
251  } else {
252  d = d2/ASTRO_JR;
253  }
254 
255  // DSM occuring at time nu2*T2
256  kep_toolbox::diff(v_out, v_beg_l, v);
257  DV[i] = kep_toolbox::norm(v_out) + std::max((2.0-d),0.0) * 1000.0;
258  }
259  // Now we return the constraints
260  c[0] = std::accumulate(T.begin(),T.end(),0.0) - m_tmax;
261  c[1] = m_dmin - d;
262 //Here the lambert solver or the lagrangian propagator went wrong
263 } catch (...) {
264  c[0] = boost::numeric::bounds<double>::highest();
265  c[1] = boost::numeric::bounds<double>::highest();
266 }
267 }
268 
270 
279 std::string mga_incipit_cstrs::pretty(const std::vector<double> &x) const {
280 
281  // We set the std output format
282  std::ostringstream s;
283  s.precision(15);
284  s << std::scientific;
285 
286  double d,ra,d2,ra2;
287 
288  double common_mu = m_seq[0]->get_mu_central_body();
289  // 1 - we 'decode' the chromosome recording the various times of flight (days) in the list T
290  std::vector<double> T(m_seq.size(),0.0);
291 
292  for (size_t i = 0; i<m_seq.size(); ++i) {
293  T[i] = x[4*i+3];
294  }
295  // 2 - We compute the epochs and ephemerides of the planetary encounters
296  std::vector<kep_toolbox::epoch> t_P(m_seq.size());
297  std::vector<kep_toolbox::array3D> r_P(m_seq.size());
298  std::vector<kep_toolbox::array3D> v_P(m_seq.size());
299  std::vector<double> DV(m_seq.size());
300  for (size_t i = 0; i<r_P.size(); ++i) {
301  t_P[i] = kep_toolbox::epoch(x[0] + std::accumulate(T.begin(),T.begin()+1+i,0.0));
302  m_seq[i]->eph(t_P[i], r_P[i], v_P[i]);
303  }
304 
305  // 3 - We start with the first leg
306  double theta = 2*boost::math::constants::pi<double>()*x[1];
307  double phi = acos(2*x[2]-1)-boost::math::constants::pi<double>() / 2;
308  kep_toolbox::array3D r = { {ASTRO_JR * 1000*cos(phi)*sin(theta), ASTRO_JR * 1000*cos(phi)*cos(theta), ASTRO_JR * 1000*sin(phi)} };
309  kep_toolbox::array3D v;
310 
311  kep_toolbox::lambert_problem l(r,r_P[0],T[0]*ASTRO_DAY2SEC,common_mu,false,false);
312  kep_toolbox::array3D v_beg_l = l.get_v1()[0];
313  kep_toolbox::array3D v_end_l = l.get_v2()[0];
314  kep_toolbox::closest_distance(d,ra,r,v_beg_l, r_P[0], v_end_l, common_mu);
315 
316  DV[0] = std::abs(kep_toolbox::norm(v_beg_l)-3400.0);
317  kep_toolbox::array3D v_out,mem_vin,mem_vout,mem_vP;
318 
319  s << "\nFirst Leg: 1000JR to " << m_seq[0]->get_name() << std::endl;
320  s << "\tDeparture: " << kep_toolbox::epoch(x[0]) << " (" << x[0] << " mjd2000) " << std::endl;
321  s << "\tDuration: " << T[0] << "days" << std::endl;
322  s << "\tInitial Velocity Increment (m/s): " << DV[0] << std::endl;
323  kep_toolbox::diff(v_out, v_end_l, v_P[0]);
324  s << "\tArrival relative velocity at " << m_seq[0]->get_name() << " (m/s): " << kep_toolbox::norm(v_out) << std::endl;
325  s << "\tClosest distance: " << d / ASTRO_JR;
326 
327  // 4 - And we proceed with each successive leg (if any)
328  for (size_t i = 1; i<m_seq.size(); ++i) {
329  // Fly-by
330  kep_toolbox::fb_prop(v_out, v_end_l, v_P[i-1], x[4*i+1] * m_seq[i-1]->get_radius(), x[4*i], m_seq[i-1]->get_mu_self());
331  // s/c propagation before the DSM
332  r = r_P[i-1];
333  v = v_out;
334  mem_vout = v_out;
335  mem_vin = v_end_l;
336  mem_vP = v_P[i-1];
337 
338  kep_toolbox::propagate_lagrangian(r,v,x[4*i+2]*T[i]*ASTRO_DAY2SEC,common_mu);
339  kep_toolbox::closest_distance(d, ra, r_P[i-1], v_out, r, v, common_mu);
340 
341  // Lambert arc to reach Earth during (1-nu2)*T2 (second segment)
342  double dt = (1-x[4*i+2])*T[i]*ASTRO_DAY2SEC;
343  kep_toolbox::lambert_problem l2(r,r_P[i],dt,common_mu,false,false);
344  v_end_l = l2.get_v2()[0];
345  v_beg_l = l2.get_v1()[0];
346  kep_toolbox::closest_distance(d2,ra2,r,v_beg_l, r_P[i], v_end_l, common_mu);
347 
348  if (d < d2)
349  {
350  d = d/ASTRO_JR;
351  ra = ra/ASTRO_JR;
352  } else {
353  d = d2/ASTRO_JR;
354  ra = ra2/ASTRO_JR;
355  }
356  // DSM occuring at time nu2*T2
357  kep_toolbox::diff(v_out, v_beg_l, v);
358  DV[i] = kep_toolbox::norm(v_out);
359  s << "\nleg no. " << i+1 << ": " << m_seq[i-1]->get_name() << " to " << m_seq[i]->get_name() << std::endl;
360  s << "\tDuration: (days)" << T[i] << std::endl;
361  s << "\tFly-by epoch: " << t_P[i-1] << " (" << t_P[i-1].mjd2000() << " mjd2000) " << std::endl;
362  s << "\tFly-by altitude (km): " << (x[4*i+1]*m_seq[i-1]->get_radius()-m_seq[i-1]->get_radius())/1000.0 << std::endl;
363  s << "\tPlanet position (m): " << r_P[i-1] << std::endl;
364  s << "\tPlanet velocity (m/s): " << mem_vP << std::endl;
365  s << "\tV in (m/s): " << mem_vin << std::endl;
366  s << "\tV out (m/s): " << mem_vout << std::endl << std::endl;
367 
368  s << "\tDSM after (days): " << x[4*i+2]*T[i] << std::endl;
369  s << "\tDSM magnitude (m/s): " << DV[i] << std::endl;
370  s << "\tClosest distance (JR): " << d << std::endl;
371  s << "\tApoapsis at closest distance (JR): " << ra << std::endl;
372  }
373 
374  s << "\nArrival at " << m_seq[m_seq.size()-1]->get_name() << std::endl;
375  kep_toolbox::diff(v_out, v_end_l, v_P[m_seq.size()-1]);
376  s << "Arrival epoch: " << t_P[m_seq.size()-1] << " (" << t_P[m_seq.size()-1].mjd2000() << " mjd2000) " << std::endl;
377  s << "Arrival Vinf (m/s): " << v_out << " - " << kep_toolbox::norm(v_out) << std::endl;
378  s << "Total mission time (days): " << std::accumulate(T.begin(),T.end(),0.0) << std::endl;
379  return s.str();
380 }
381 std::string mga_incipit_cstrs::get_name() const
382 {
383  return "MGA-INCIPIT (CAPTURE AT JUPITER) - Constrained version";
384 }
385 
386 
388 
391 const std::vector<std::vector<double> >& mga_incipit_cstrs::get_tof() const {
392  return m_tof;
393 }
394 
396 
401 void mga_incipit_cstrs::set_tof(const std::vector<std::vector<double> >& tof) {
402  if (tof.size() != (m_seq.size())) {
403  pagmo_throw(value_error,"The time-of-flight vector (tof) has the wrong length");
404  }
405  m_tof = tof;
406  for (size_t i=0; i< m_seq.size(); ++i) {
407  set_bounds(3+4*i,tof[i][0],tof[i][1]);
408  }
409 }
410 
412 
415 std::vector<kep_toolbox::planet::planet_ptr> mga_incipit_cstrs::get_sequence() const {
416  return m_seq;
417 }
418 
420 
425 {
426  std::ostringstream oss;
427  oss << "\n\tSequence: ";
428  for (size_t i = 0; i<m_seq.size() ;++i) {
429  oss << m_seq[i]->get_name() << " ";
430  }
431  oss << "\n\tTime of flights?: ";
432  for (size_t i=0; i<m_seq.size(); ++i) {
433  oss << m_tof[i]<<' ';
434  }
435  oss << "\n\tMaximum time of flight: " << m_tmax;
436  oss << "\n\tMinimum distance from Jupiter: " << m_dmin;
437 
438  return oss.str();
439 }
440 
441 }} //namespaces
442 
443 BOOST_CLASS_EXPORT_IMPLEMENT(pagmo::problem::mga_incipit_cstrs)
Root PaGMO namespace.
void objfun_impl(fitness_vector &, const decision_vector &) const
Implementation of the objective function.
boost::shared_ptr< base > base_ptr
Alias for shared pointer to base problem.
Definition: problem/base.h:62
std::string get_name() const
Get problem's name.
std::vector< double > decision_vector
Decision vector type.
Definition: types.h:40
void set_tof(const std::vector< std::vector< double > > &)
Sets the times of flight.
Base problem class.
Definition: problem/base.h:148
mga_incipit_cstrs(const std::vector< kep_toolbox::planet::planet_ptr >=construct_default_sequence(), const kep_toolbox::epoch t0_l=kep_toolbox::epoch(7305.0), const kep_toolbox::epoch t0_u=kep_toolbox::epoch(11323.0), const std::vector< std::vector< double > > tof=construct_default_tofs(), double Tmax=365.25, double Dmin=0.2)
Constructor.
base_ptr clone() const
Clone method.
std::string human_readable_extra() const
Extra human readable info for the problem.
void compute_constraints_impl(constraint_vector &, const decision_vector &) const
Implementation of the objective function.
The beginning of the GTOC6 Jupiter Capture Trajectory.
std::vector< double > fitness_vector
Fitness vector type.
Definition: types.h:42
std::vector< double > constraint_vector
Constraint vector type.
Definition: types.h:44
const decision_vector & get_ub() const
Upper bounds getter.
const std::vector< std::vector< double > > & get_tof() const
Gets the times of flight.
const decision_vector & get_lb() const
Lower bounds getter.
std::vector< kep_toolbox::planet::planet_ptr > get_sequence() const
Gets the planetary sequence defining the interplanetary mga-1dsm mission.
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
std::string pretty(const std::vector< double > &x) const
Outputs a stream with the trajectory data.