PaGMO  1.1.5
tandem.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 <keplerian_toolbox/epoch.h>
27 
28 #include "tandem.h"
29 #include "../AstroToolbox/mga_dsm.h"
30 #include "../AstroToolbox/misc4Tandem.h"
31 
32 namespace pagmo { namespace problem {
33 
35 const int tandem::Data[24][5] = {
36  {3,2,2,2,6},
37  {3,2,2,3,6},
38  {3,2,2,4,6},
39  {3,2,2,5,6},
40  {3,2,3,2,6},
41  {3,2,3,3,6},
42  {3,2,3,4,6},
43  {3,2,3,5,6},
44  {3,2,4,2,6},
45  {3,2,4,3,6},
46  {3,2,4,4,6},
47  {3,2,4,5,6},
48  {3,3,2,2,6},
49  {3,3,2,3,6},
50  {3,3,2,4,6},
51  {3,3,2,5,6},
52  {3,3,3,2,6},
53  {3,3,3,3,6},
54  {3,3,3,4,6},
55  {3,3,3,5,6},
56  {3,3,4,2,6},
57  {3,3,4,3,6},
58  {3,3,4,4,6},
59  {3,3,4,5,6}
60 };
61 
63 const int tandem::sequence[5] = {1,1,1,1,1};
64 
66 
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)
72 {
73  if (probid < 1 || probid > 24) {
74  pagmo_throw(value_error,"probid needs to be an integer in [1,24]");
75  }
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)");
78  }
79  for (int i =0; i<5; i++) {
80  problem.sequence[i] = Data[probid-1][i];
81  }
82 
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};
85 
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};
88 
89  if (tof_==-1) set_bounds(lbunc,lbunc+18,ubunc,ubunc+18);
90  else set_bounds(lbcon,lbcon+18,ubcon,ubcon+18);
91 }
92 
95 {
96  return base_ptr(new tandem(*this));
97 }
98 
101 {
102  if (tof!=-1) { //constrained problem
103  //Here we copy the chromosome into a new vector and we transform its time percentages into days
104  copy_of_x = x;
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]);
110  } else { //unconstrained problem
111  MGA_DSM(x, problem, f[0]);
112  }
113  //evaluating the mass from the dvs
114  double rE[3];
115  double vE[3];
116  Planet_Ephemerides_Analytical (x[0],3,rE,vE);
117  double VINFE = x[1];
118  double udir = x[2];
119  double vdir = x[3];
120  double vtemp[3];
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];
124  double iP1[3];
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;
129  double zP1[3];
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;
134  double jP1[3];
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; //See Picking a Point on a Sphere
139  double phi=acos(2*vdir-1)-M_PI/2; //In this way: -pi/2<phi<pi/2 so phi can be used as out-of-plane rotation
140  double vinf[3];
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]);
144  //We rotate it to the equatorial plane
145  ecl2equ(vinf,vinf);
146  //And we find the declination in degrees
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;
150 
151  //double m_initial = SoyuzFregat(VINFE,declination);
152  double m_initial = Atlas501(VINFE,declination);
153 
154  //We evaluate the final mass
155  double Isp = 312;
156  double g0 = 9.80665;
157  double sumDVvec=0;
158 
159  for(unsigned int i=1;i<=5;i++) {
160  sumDVvec=sumDVvec+problem.DV[i];
161  }
162  double m_final;
163  sumDVvec=sumDVvec+0.165; //losses for 3 swgbys + insertion
164  m_final = m_initial * exp(-sumDVvec/Isp/g0*1000);
165  f[0] = -log(m_final);
166 }
167 
169 
177 std::string tandem::pretty(const std::vector<double> &x) const
178 {
179  double obj=0;
180  std::vector<double> printablex;
181  if (tof!=-1) { //constrained problem
182  //Here we copy the chromosome into a new vector and we transform its time percentages into days
183  copy_of_x = x;
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;
190  } else { //unconstrained problem
191  MGA_DSM(x, problem, obj);
192  printablex=x;
193  }
194  std::ostringstream s;
195  s.precision(15);
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];
203  }
204  s << '\n';
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] << ' ';
210  }
211  s << '\n';
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];
216  }
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';
220  }
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';
225  }
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';
228  }
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';
231  }
232  s << "Final DV:\t\t\t" << problem.DV.back() << '\n';
233  return s.str();
234 }
235 
237 
244 void tandem::set_sparsity(int &lenG, std::vector<int> &iGfun, std::vector<int> &jGvar) const
245 {
246  lenG=18;
247  iGfun.resize(18);
248  jGvar.resize(18);
249  for (int i = 0; i<lenG; ++i)
250  {
251  iGfun[i] = 0;
252  jGvar[i] = i;
253  }
254 }
255 
256 std::string tandem::get_name() const
257 {
258  return "TandEM";
259 }
260 
261 }} //namespaces
262 
263 BOOST_CLASS_EXPORT_IMPLEMENT(pagmo::problem::tandem)
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
Base problem class.
Definition: problem/base.h:148
base_ptr clone() const
Clone method.
Definition: tandem.cpp:94
std::string get_name() const
Get problem's name.
Definition: tandem.cpp:256
std::vector< double > fitness_vector
Fitness vector type.
Definition: types.h:42
void set_sparsity(int &, std::vector< int > &, std::vector< int > &) const
Implementation of the sparsity structure.
Definition: tandem.cpp:244
tandem(const int problemid=6, const double tof_=-1)
Constructor.
Definition: tandem.cpp:71
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.
Definition: tandem.cpp:100
std::string pretty(const std::vector< double > &x) const
Outputs a stream of the trajectory data.
Definition: tandem.cpp:177
TandEM problem.
Definition: tandem.h:80