26 #include <keplerian_toolbox/epoch.h>
29 #include "../AstroToolbox/mga_dsm.h"
30 #include "../AstroToolbox/misc4Tandem.h"
32 namespace pagmo {
namespace problem {
35 const int tandem::Data[24][5] = {
63 const int tandem::sequence[5] = {1,1,1,1,1};
71 tandem::tandem(
const int probid,
const double tof_):
base(18), problem(orbit_insertion,sequence,5,0,0,0,0.98531407996358,80330.0), tof(tof_),copy_of_x(18)
73 if (probid < 1 || probid > 24) {
74 pagmo_throw(value_error,
"probid needs to be an integer in [1,24]");
76 if (tof_ > 20 && tof_ <5 && tof_!=-1) {
77 pagmo_throw(value_error,
"time of flight constraint needs to be a number in [5,20] (years)");
79 for (
int i =0; i<5; i++) {
80 problem.sequence[i] = Data[probid-1][i];
83 const double lbunc[18] = {5475, 2.50001, 0, 0, 20 , 20 , 20 , 20 , 0.01, 0.01, 0.01, 0.01, 1.05, 1.05, 1.05, -M_PI, -M_PI, -M_PI};
84 const double ubunc[18] = {9132, 4.9, 1, 1, 2500 , 2500, 2500, 2500, 0.99, 0.99, 0.99, 0.99, 10, 10, 10, M_PI, M_PI, M_PI};
86 const double lbcon[18] = {5475, 2.50001, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 1.05, 1.05, 1.05, -M_PI, -M_PI, -M_PI};
87 const double ubcon[18] = {9132, 4.9, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 0.99, 10, 10, 10, M_PI, M_PI, M_PI};
89 if (tof_==-1)
set_bounds(lbunc,lbunc+18,ubunc,ubunc+18);
90 else set_bounds(lbcon,lbcon+18,ubcon,ubcon+18);
105 copy_of_x[4] = x[4]*365.25*tof;
106 copy_of_x[5] = x[5]*(365.25*tof-copy_of_x[4]);
107 copy_of_x[6] = x[6]*(365.25*tof-copy_of_x[4]-copy_of_x[5]);
108 copy_of_x[7] = x[7]*(365.25*tof-copy_of_x[4]-copy_of_x[5]-copy_of_x[6]);
109 MGA_DSM(copy_of_x, problem, f[0]);
111 MGA_DSM(x, problem, f[0]);
116 Planet_Ephemerides_Analytical (x[0],3,rE,vE);
121 vtemp[0]= rE[1]*vE[2]-rE[2]*vE[1];
122 vtemp[1]= rE[2]*vE[0]-rE[0]*vE[2];
123 vtemp[2]= rE[0]*vE[1]-rE[1]*vE[0];
125 double normvE=sqrt(vE[0]*vE[0]+vE[1]*vE[1]+vE[2]*vE[2]);
126 iP1[0]= vE[0]/normvE;
127 iP1[1]= vE[1]/normvE;
128 iP1[2]= vE[2]/normvE;
130 double normvtemp=sqrt(vtemp[0]*vtemp[0]+vtemp[1]*vtemp[1]+vtemp[2]*vtemp[2]);
131 zP1[0]= vtemp[0]/normvtemp;
132 zP1[1]= vtemp[1]/normvtemp;
133 zP1[2]= vtemp[2]/normvtemp;
135 jP1[0]= zP1[1]*iP1[2]-zP1[2]*iP1[1];
136 jP1[1]= zP1[2]*iP1[0]-zP1[0]*iP1[2];
137 jP1[2]= zP1[0]*iP1[1]-zP1[1]*iP1[0];
138 double theta=2*M_PI*udir;
139 double phi=acos(2*vdir-1)-M_PI/2;
141 vinf[0]=VINFE*(cos(theta)*cos(phi)*iP1[0]+sin(theta)*cos(phi)*jP1[0]+sin(phi)*zP1[0]);
142 vinf[1]=VINFE*(cos(theta)*cos(phi)*iP1[1]+sin(theta)*cos(phi)*jP1[1]+sin(phi)*zP1[1]);
143 vinf[2]=VINFE*(cos(theta)*cos(phi)*iP1[2]+sin(theta)*cos(phi)*jP1[2]+sin(phi)*zP1[2]);
147 double normvinf=sqrt(vinf[0]*vinf[0]+vinf[1]*vinf[1]+vinf[2]*vinf[2]);
148 double sindelta = vinf[2] / normvinf;
149 double declination = asin(sindelta)/M_PI*180;
152 double m_initial = Atlas501(VINFE,declination);
159 for(
unsigned int i=1;i<=5;i++) {
160 sumDVvec=sumDVvec+problem.DV[i];
163 sumDVvec=sumDVvec+0.165;
164 m_final = m_initial * exp(-sumDVvec/Isp/g0*1000);
165 f[0] = -log(m_final);
180 std::vector<double> printablex;
184 copy_of_x[4] = x[4]*365.25*tof;
185 copy_of_x[5] = x[5]*(365.25*tof-copy_of_x[4]);
186 copy_of_x[6] = x[6]*(365.25*tof-copy_of_x[4]-copy_of_x[5]);
187 copy_of_x[7] = x[7]*(365.25*tof-copy_of_x[4]-copy_of_x[5]-copy_of_x[6]);
188 MGA_DSM(copy_of_x, problem, obj);
189 printablex=copy_of_x;
191 MGA_DSM(x, problem, obj);
194 std::ostringstream s;
196 s << std::scientific;
197 const size_t seq_size = (x.size() + 2) / 4;
198 pagmo_assert((x.size() + 2) % 4 == 0 && seq_size >= 2);
199 pagmo_assert(problem.sequence.size() == seq_size);
200 s <<
"Flyby sequence:\t\t\t";
201 for (
size_t i = 0; i < seq_size; ++i) {
202 s << problem.sequence[i];
205 s <<
"Departure epoch (mjd2000):\t" << printablex[0] <<
'\n';
206 s <<
"Departure epoch:\t\t" << ::kep_toolbox::epoch(printablex[0],::kep_toolbox::epoch::MJD2000) <<
'\n';
207 s <<
"Vinf polar components:\t\t";
208 for (
size_t i = 0; i < 3; ++i) {
209 s << printablex[i + 1] <<
' ';
212 double totaltime = 0;
213 for (
size_t i = 0; i < seq_size - 1; ++i) {
214 s <<
"Leg time of flight:\t\t" << printablex[i + 4] <<
'\n';
215 totaltime += printablex[i + 4];
217 s <<
"Total time of flight:\t\t" << totaltime <<
'\n';
218 for (
size_t i = 0; i < seq_size - 2; ++i) {
219 s <<
"Flyby radius:\t\t\t" << printablex[i + 2 * (seq_size + 1)] <<
'\n';
221 totaltime=printablex[0];
222 for (
size_t i = 0; i < seq_size - 2; ++i) {
223 totaltime += printablex[i + 4];
224 s <<
"Flyby date:\t\t\t" << ::kep_toolbox::epoch(totaltime,::kep_toolbox::epoch::MJD2000) <<
'\n';
226 for (
size_t i = 0; i < seq_size - 2; ++i) {
227 s <<
"Vinf at flyby:\t\t\t" << std::sqrt(problem.vrelin_vec[i]) <<
'\n';
229 for (
size_t i = 0; i < seq_size - 1; ++i) {
230 s <<
"dsm" << i+1 <<
":\t\t\t\t" << problem.DV[i+1] <<
'\n';
232 s <<
"Final DV:\t\t\t" << problem.DV.back() <<
'\n';
249 for (
int i = 0; i<lenG; ++i)
boost::shared_ptr< base > base_ptr
Alias for shared pointer to base problem.
std::vector< double > decision_vector
Decision vector type.
base_ptr clone() const
Clone method.
std::string get_name() const
Get problem's name.
std::vector< double > fitness_vector
Fitness vector type.
void set_sparsity(int &, std::vector< int > &, std::vector< int > &) const
Implementation of the sparsity structure.
tandem(const int problemid=6, const double tof_=-1)
Constructor.
void set_bounds(const decision_vector &, const decision_vector &)
Bounds setter from pagmo::decision_vector.
void objfun_impl(fitness_vector &, const decision_vector &) const
Implementation of the objective function.
std::string pretty(const std::vector< double > &x) const
Outputs a stream of the trajectory data.