26 #include <boost/math/constants/constants.hpp>
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>
35 #include "mga_incipit_cstrs.h"
38 namespace pagmo {
namespace problem {
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,
62 double Dmin) :
base(4*seq.size(),0,1,2,2,1E-3), m_tof(tof), m_tmax(Tmax), m_dmin(Dmin)
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();
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");
73 if (tof.size() != seq.size()) {
74 pagmo_throw(value_error,
"The time-of-flight vector (tof) has the wrong length");
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)");
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());
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];
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];
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();
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()),
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());
133 double common_mu = m_seq[0]->get_mu_central_body();
135 std::vector<double> T(m_seq.size(),0.0);
137 for (
size_t i = 0; i<m_seq.size(); ++i) {
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]);
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;
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];
160 DV[0] = std::abs(kep_toolbox::norm(v_beg_l)-3400.0);
163 kep_toolbox::array3D v_out;
164 for (
size_t i = 1; i<m_seq.size(); ++i) {
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());
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);
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);
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;
191 f[0] = std::accumulate(DV.begin(),DV.end(),0.0);
194 f[0] = boost::numeric::bounds<double>::highest();
202 double common_mu = m_seq[0]->get_mu_central_body();
204 std::vector<double> T(m_seq.size(),0.0);
206 for (
size_t i = 0; i<m_seq.size(); ++i) {
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]);
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;
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];
229 DV[0] = std::abs(kep_toolbox::norm(v_beg_l)-3400.0);
232 kep_toolbox::array3D v_out;
233 for (
size_t i = 1; i<m_seq.size(); ++i) {
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());
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);
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);
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;
260 c[0] = std::accumulate(T.begin(),T.end(),0.0) - m_tmax;
264 c[0] = boost::numeric::bounds<double>::highest();
265 c[1] = boost::numeric::bounds<double>::highest();
282 std::ostringstream s;
284 s << std::scientific;
288 double common_mu = m_seq[0]->get_mu_central_body();
290 std::vector<double> T(m_seq.size(),0.0);
292 for (
size_t i = 0; i<m_seq.size(); ++i) {
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]);
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;
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);
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;
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;
328 for (
size_t i = 1; i<m_seq.size(); ++i) {
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());
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);
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);
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;
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;
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;
383 return "MGA-INCIPIT (CAPTURE AT JUPITER) - Constrained version";
402 if (tof.size() != (m_seq.size())) {
403 pagmo_throw(value_error,
"The time-of-flight vector (tof) has the wrong length");
406 for (
size_t i=0; i< m_seq.size(); ++i) {
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() <<
" ";
431 oss <<
"\n\tTime of flights?: ";
432 for (
size_t i=0; i<m_seq.size(); ++i) {
433 oss << m_tof[i]<<
' ';
435 oss <<
"\n\tMaximum time of flight: " << m_tmax;
436 oss <<
"\n\tMinimum distance from Jupiter: " << m_dmin;
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.
std::string get_name() const
Get problem's name.
std::vector< double > decision_vector
Decision vector type.
void set_tof(const std::vector< std::vector< double > > &)
Sets the times of flight.
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.
std::vector< double > constraint_vector
Constraint vector type.
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.
std::string pretty(const std::vector< double > &x) const
Outputs a stream with the trajectory data.