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.h"
38 namespace pagmo {
namespace problem {
55 const kep_toolbox::epoch t0_l,
const kep_toolbox::epoch t0_u,
56 const std::vector<std::vector<double> > tof
58 base(4*seq.size()), m_tof(tof)
61 std::vector<double> mus(seq.size());
62 for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i< seq.size(); ++i) {
63 mus[i] = seq[i]->get_mu_central_body();
65 if ((
unsigned int)std::count(mus.begin(), mus.end(), mus[0]) != mus.size()) {
66 pagmo_throw(value_error,
"The planets do not all have the same mu_central_body");
69 if (tof.size() != seq.size()) {
70 pagmo_throw(value_error,
"The time-of-flight vector (tof) has the wrong length");
72 for (
size_t i = 0; i < tof.size(); ++i) {
73 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)");
77 for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i < seq.size(); ++i) {
78 m_seq.push_back(seq[i]->
clone());
86 lb[0] = t0_l.mjd2000(); ub[0] = t0_u.mjd2000();
87 lb[1] = 0; lb[2] = 0; ub[1] = 1; ub[2] = 1;
88 lb[3] = m_tof[0][0]; ub[3] = m_tof[0][1];
91 for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 1; i < m_tof.size(); ++i) {
92 lb[4*i] = - 2 * boost::math::constants::pi<double>(); ub[4*i] = 2 * boost::math::constants::pi<double>();
93 lb[4*i+1] = 1.1; ub[4*i+1] = 30;
94 lb[4*i+2] = 1e-5; ub[4*i+2] = 1-1e-5;
95 lb[4*i+3] = m_tof[i][0]; ub[4*i+3] = m_tof[i][1];
99 for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i < m_tof.size()-1; ++i) {
100 lb[4*i+5] = m_seq[i]->get_safe_radius() / m_seq[i]->get_radius();
101 ub[4*i+5] = (m_seq[i]->get_radius() + 2000000) / m_seq[i]->get_radius();
109 for (std::vector<kep_toolbox::planet::planet_ptr>::size_type i = 0; i < p.m_seq.size();++i) {
110 m_seq.push_back(p.m_seq[i]->clone());
125 double common_mu = m_seq[0]->get_mu_central_body();
127 std::vector<double> T(m_seq.size(),0.0);
129 for (
size_t i = 0; i<m_seq.size(); ++i) {
133 std::vector<kep_toolbox::epoch> t_P(m_seq.size());
134 std::vector<kep_toolbox::array3D> r_P(m_seq.size());
135 std::vector<kep_toolbox::array3D> v_P(m_seq.size());
136 std::vector<double> DV(m_seq.size());
137 for (
size_t i = 0; i<r_P.size(); ++i) {
138 t_P[i] = kep_toolbox::epoch(x[0] + std::accumulate(T.begin(),T.begin()+1+i,0.0));
139 m_seq[i]->eph(t_P[i], r_P[i], v_P[i]);
143 double theta = 2*boost::math::constants::pi<double>()*x[1];
144 double phi = acos(2*x[2]-1)-boost::math::constants::pi<double>() / 2;
146 kep_toolbox::array3D r = { {ASTRO_JR*1000*cos(phi)*sin(theta), ASTRO_JR*1000*cos(phi)*cos(theta), ASTRO_JR*1000*sin(phi)} };
147 kep_toolbox::array3D v;
148 kep_toolbox::lambert_problem l(r,r_P[0],T[0]*ASTRO_DAY2SEC,common_mu,
false,
false);
149 kep_toolbox::array3D v_beg_l = l.get_v1()[0];
150 kep_toolbox::array3D v_end_l = l.get_v2()[0];
152 DV[0] = std::abs(kep_toolbox::norm(v_beg_l)-3400.0);
155 kep_toolbox::array3D v_out;
156 for (
size_t i = 1; i<m_seq.size(); ++i) {
158 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());
162 kep_toolbox::propagate_lagrangian(r,v,x[4*i+2]*T[i]*ASTRO_DAY2SEC,common_mu);
163 kep_toolbox::closest_distance(d, ra, r_P[i-1], v_out, r, v, common_mu);
166 double dt = (1-x[4*i+2])*T[i]*ASTRO_DAY2SEC;
167 kep_toolbox::lambert_problem l2(r,r_P[i],dt,common_mu,
false,
false);
168 v_end_l = l2.get_v2()[0];
169 v_beg_l = l2.get_v1()[0];
170 kep_toolbox::closest_distance(d2,ra2,r,v_beg_l, r_P[i], v_end_l, common_mu);
179 kep_toolbox::diff(v_out, v_beg_l, v);
180 DV[i] = kep_toolbox::norm(v_out) + std::max((2.0-d),0.0) * 1000.0;
183 f[0] = std::accumulate(DV.begin(),DV.end(),0.0);
186 f[0] = boost::numeric::bounds<double>::highest();
203 std::ostringstream s;
205 s << std::scientific;
209 double common_mu = m_seq[0]->get_mu_central_body();
211 std::vector<double> T(m_seq.size(),0.0);
213 for (
size_t i = 0; i<m_seq.size(); ++i) {
217 std::vector<kep_toolbox::epoch> t_P(m_seq.size());
218 std::vector<kep_toolbox::array3D> r_P(m_seq.size());
219 std::vector<kep_toolbox::array3D> v_P(m_seq.size());
220 std::vector<double> DV(m_seq.size());
221 for (
size_t i = 0; i<r_P.size(); ++i) {
222 t_P[i] = kep_toolbox::epoch(x[0] + std::accumulate(T.begin(),T.begin()+1+i,0.0));
223 m_seq[i]->eph(t_P[i], r_P[i], v_P[i]);
227 double theta = 2*boost::math::constants::pi<double>()*x[1];
228 double phi = acos(2*x[2]-1)-boost::math::constants::pi<double>() / 2;
229 kep_toolbox::array3D r = { {ASTRO_JR * 1000*cos(phi)*sin(theta), ASTRO_JR * 1000*cos(phi)*cos(theta), ASTRO_JR * 1000*sin(phi)} };
230 kep_toolbox::array3D v;
232 kep_toolbox::lambert_problem l(r,r_P[0],T[0]*ASTRO_DAY2SEC,common_mu,
false,
false);
233 kep_toolbox::array3D v_beg_l = l.get_v1()[0];
234 kep_toolbox::array3D v_end_l = l.get_v2()[0];
235 kep_toolbox::closest_distance(d,ra,r,v_beg_l, r_P[0], v_end_l, common_mu);
237 DV[0] = std::abs(kep_toolbox::norm(v_beg_l)-3400.0);
238 kep_toolbox::array3D v_out,mem_vin,mem_vout,mem_vP;
240 s <<
"\nFirst Leg: 1000JR to " << m_seq[0]->get_name() << std::endl;
241 s <<
"\tDeparture: " << kep_toolbox::epoch(x[0]) <<
" (" << x[0] <<
" mjd2000) " << std::endl;
242 s <<
"\tDuration: " << T[0] <<
"days" << std::endl;
243 s <<
"\tInitial Velocity Increment (m/s): " << DV[0] << std::endl;
244 kep_toolbox::diff(v_out, v_end_l, v_P[0]);
245 s <<
"\tArrival relative velocity at " << m_seq[0]->get_name() <<
" (m/s): " << kep_toolbox::norm(v_out) << std::endl;
246 s <<
"\tClosest distance: " << d / ASTRO_JR;
249 for (
size_t i = 1; i<m_seq.size(); ++i) {
251 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());
259 kep_toolbox::propagate_lagrangian(r,v,x[4*i+2]*T[i]*ASTRO_DAY2SEC,common_mu);
260 kep_toolbox::closest_distance(d, ra, r_P[i-1], v_out, r, v, common_mu);
263 double dt = (1-x[4*i+2])*T[i]*ASTRO_DAY2SEC;
264 kep_toolbox::lambert_problem l2(r,r_P[i],dt,common_mu,
false,
false);
265 v_end_l = l2.get_v2()[0];
266 v_beg_l = l2.get_v1()[0];
267 kep_toolbox::closest_distance(d2,ra2,r,v_beg_l, r_P[i], v_end_l, common_mu);
278 kep_toolbox::diff(v_out, v_beg_l, v);
279 DV[i] = kep_toolbox::norm(v_out);
280 s <<
"\nleg no. " << i+1 <<
": " << m_seq[i-1]->get_name() <<
" to " << m_seq[i]->get_name() << std::endl;
281 s <<
"\tDuration: (days)" << T[i] << std::endl;
282 s <<
"\tFly-by epoch: " << t_P[i-1] <<
" (" << t_P[i-1].mjd2000() <<
" mjd2000) " << std::endl;
283 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;
284 s <<
"\tPlanet position (m): " << r_P[i-1] << std::endl;
285 s <<
"\tPlanet velocity (m/s): " << mem_vP << std::endl;
286 s <<
"\tV in (m/s): " << mem_vin << std::endl;
287 s <<
"\tV out (m/s): " << mem_vout << std::endl << std::endl;
289 s <<
"\tDSM after (days): " << x[4*i+2]*T[i] << std::endl;
290 s <<
"\tDSM magnitude (m/s): " << DV[i] << std::endl;
291 s <<
"\tClosest distance (JR): " << d << std::endl;
292 s <<
"\tApoapsis at closest distance (JR): " << ra << std::endl;
295 s <<
"\nArrival at " << m_seq[m_seq.size()-1]->get_name() << std::endl;
296 kep_toolbox::diff(v_out, v_end_l, v_P[m_seq.size()-1]);
297 s <<
"Arrival epoch: " << t_P[m_seq.size()-1] <<
" (" << t_P[m_seq.size()-1].mjd2000() <<
" mjd2000) " << std::endl;
298 s <<
"Arrival Vinf (m/s): " << v_out <<
" - " << kep_toolbox::norm(v_out) << std::endl;
299 s <<
"Total mission time (days): " << std::accumulate(T.begin(),T.end(),0.0) << std::endl;
304 return "MGA-INCIPIT (CAPTURE AT JUPITER)";
323 if (tof.size() != (m_seq.size())) {
324 pagmo_throw(value_error,
"The time-of-flight vector (tof) has the wrong length");
327 for (
size_t i=0; i< m_seq.size(); ++i) {
347 std::ostringstream oss;
348 oss <<
"\n\tSequence: ";
349 for (
size_t i = 0; i<m_seq.size() ;++i) {
350 oss << m_seq[i]->get_name() <<
" ";
352 oss <<
"\n\tTime of flights?: ";
353 for (
size_t i=0; i<m_seq.size(); ++i) {
354 oss << m_tof[i]<<
' ';
boost::shared_ptr< base > base_ptr
Alias for shared pointer to base problem.
std::vector< double > decision_vector
Decision vector type.
std::string get_name() const
Get problem's name.
std::string pretty(const std::vector< double > &x) const
Outputs a stream with the trajectory data.
void set_tof(const std::vector< std::vector< double > > &)
Sets the times of flight.
void objfun_impl(fitness_vector &, const decision_vector &) const
Implementation of the objective function.
const std::vector< std::vector< double > > & get_tof() const
Gets the times of flight.
std::vector< kep_toolbox::planet::planet_ptr > get_sequence() const
Gets the planetary sequence defining the interplanetary mga-1dsm mission.
std::vector< double > fitness_vector
Fitness vector type.
const decision_vector & get_ub() const
Upper bounds getter.
base_ptr clone() const
Clone method.
const decision_vector & get_lb() const
Lower bounds getter.
void set_bounds(const decision_vector &, const decision_vector &)
Bounds setter from pagmo::decision_vector.
The beginning of the GTOC6 Jupiter Capture Trajectory.
std::string human_readable_extra() const
Extra human readable info for the problem.
decision_vector::size_type size_type
Problem's size type: the same as pagmo::decision_vector's size type.
mga_incipit(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())
Constructor.