Python PyKEP類代碼示例

def lambert_leg(P1, P2, i, j, t1, t2, tof, vrel=None, dv_launch=0.):
    """Compute a lambert leg from planet to planet.
    Arguments:
    p1 -- starting planet (str or PyKEP.planet object)
    p2 -- final planet (str or PyKEP.planet object)
    t0 -- start time of leg in MJD2000
    tof -- time of flight in days
    
    Keyword arguments:
    vrel -- caresian coordinates of the relative velocity before the flyby at p1
    dv_launch -- dv discounted at lunch (i.e. if vrel is None)
    rendezvous -- add final dv
    Returns:
    dV, vrel_out, where vrel_out is the relative velocity at the end of the leg at p2
    """

    ast1 = ASTEROIDS[P1]
    ast2 = ASTEROIDS[P2]

    r1 = state_asteroids.EPH[i][t1][0]
    v1 = state_asteroids.EPH[i][t1][1]
    r2 = state_asteroids.EPH[j][t2][0]
    v2 = state_asteroids.EPH[j][t2][1]

    lambert = kep.lambert_problem(r1, r2, tof * kep.DAY2SEC, ast1.mu_central_body, False, 0)

    vrel_in = tuple(map(lambda x, y: x - y, lambert.get_v1()[0], v1))
    vrel_out = tuple(map(lambda x, y: x - y, lambert.get_v2()[0], v2))

    dv_lambert = np.linalg.norm(vrel_out) + np.linalg.norm(vrel_in)

    a, _, _, dv_damon = kep.damon(vrel_in, vrel_out, tof*kep.DAY2SEC)
    m_star = kep.max_start_mass(np.linalg.norm(a), dv_damon, T_max, Isp)
    
    return dv_lambert, dv_damon, m_star

    def __init__(
            self,
            mass=1000,
            Tmax=0.05,
            Isp=2500,
            Vinf_0=3,
            Vinf_f=1e-12,
            nseg=10,
            departure=None,
            target=None,
            optimise4mass=False):
        """
        Constructs a low-thrust transfer between a departure planet and a target planet (Constrained Continuous Single-Objective)

        NOTE: An impulsive transcription is used to transform into an NLP the Optimal Control Problem

        USAGE: problem.py_pl2pl(self,mass=1000,Tmax=0.05,Isp=2500,Vinf_0=3,Vinf_f=0,nseg=10,departure = PyKEP.planet_ss('earth'), target = PyKEP.planet_ss('mars'))

        * mass: spacecraft mass at departure [kg]
        * Tmax: maximum allowed thrust [N]
        * Isp: spacecarft engine specific impulse [Isp]
        * Vinf_0: allowed maximum starting velocity [km/s]
        * Vinf_f: allowed maximum arrival velocity [km/s]
                  (if negative it is interpreted as a minimum arrival velocity)
        * nseg: number of segments used for the impulsive transcription
        * departure: departure planet (a PyKEP planet)
        * target: arrival planet (a PyKEP planet)
        """
        try:
            import PyKEP
        except ImportError:
            raise ImportError(
                "Error while trying 'import PyKEP': is PyKEP installed?")
        if departure is None:
            departure = PyKEP.planet_ss('earth')
        if target is None:
            target = mars = PyKEP.planet_ss('mars')
        super(py_pl2pl, self).__init__(
            9 + nseg * 3, 0, 1, 9 + nseg, nseg + 2, 1e-5)
        self.__departure = departure
        self.__target = target
        self.__sc = PyKEP.sims_flanagan.spacecraft(mass, Tmax, Isp)
        self.__Vinf_0 = Vinf_0 * 1000
        self.__Vinf_f = Vinf_f * 1000
        self.__leg = PyKEP.sims_flanagan.leg()
        self.__leg.set_mu(departure.mu_central_body)
        self.__leg.set_spacecraft(self.__sc)
        self.__nseg = nseg
        self.set_bounds(
            [0, 10, self.__sc.mass / 10, -abs(self.__Vinf_0),
             -abs(self.__Vinf_0), -abs(self.__Vinf_0), -abs(self.__Vinf_f),
             -abs(self.__Vinf_f), -abs(self.__Vinf_f)] + [-1] * 3 * nseg,
            [3000, 1500, self.__sc.mass, abs(self.__Vinf_0),
             abs(self.__Vinf_0), abs(self.__Vinf_0), abs(self.__Vinf_f), 
             abs(self.__Vinf_f), abs(self.__Vinf_f)] + [1] * 3 * nseg)
        self.__optimise4mass = optimise4mass

    def _compute_constraints_impl(self, x):
        import PyKEP
        start = PyKEP.epoch(x[0])
        end = PyKEP.epoch(x[0] + x[1])

        # Computing starting spaceraft state
        r, v = self.__departure.eph(start)
        v_list = list(v)
        v_list[0] += x[3]
        v_list[1] += x[4]
        v_list[2] += x[5]
        x0 = PyKEP.sims_flanagan.sc_state(r, v_list, self.__sc.mass)

        # Computing ending spaceraft state
        r, v = self.__target.eph(end)
        v_list = list(v)
        v_list[0] += x[6]
        v_list[1] += x[7]
        v_list[2] += x[8]
        xe = PyKEP.sims_flanagan.sc_state(r, v_list, x[2])

        # Building the SF leg
        self.__leg.set(start, x0, x[-3 * self.__nseg:], end, xe)

        # Computing Vinf constraints (careful here, the weights do count). In case of a larger than constarint
        # a factor of 100 has been added
        if (self.__Vinf_0 >= 0):
            v_inf_con_0 = (x[3] * x[3] + x[4] * x[4] + x[5] * x[5] - self.__Vinf_0 *
                           self.__Vinf_0) / (PyKEP.EARTH_VELOCITY * PyKEP.
                                             EARTH_VELOCITY)
        else:
            v_inf_con_0 = -100 * (x[3] * x[3] + x[4] * x[4] + x[5] * x[5] - self.
                                  __Vinf_0 * self.__Vinf_0) / (PyKEP.
                                                               EARTH_VELOCITY * PyKEP.EARTH_VELOCITY)
        if (self.__Vinf_f >= 0):
            v_inf_con_f = (x[6] * x[6] + x[7] * x[7] + x[8] * x[8] - self.__Vinf_f *
                           self.__Vinf_f) / (PyKEP.EARTH_VELOCITY * PyKEP.
                                             EARTH_VELOCITY)
        else:
            v_inf_con_f = -100 * (x[6] * x[6] + x[7] * x[7] + x[8] * x[8] - self.
                                  __Vinf_f * self.__Vinf_f) / (PyKEP.
                                                               EARTH_VELOCITY * PyKEP.EARTH_VELOCITY)

        # Setting all constraints
        retval = list(self.__leg.mismatch_constraints(
        ) + self.__leg.throttles_constraints()) + [v_inf_con_0] + [v_inf_con_f]
        retval[0] /= PyKEP.AU
        retval[1] /= PyKEP.AU
        retval[2] /= PyKEP.AU
        retval[3] /= PyKEP.EARTH_VELOCITY
        retval[4] /= PyKEP.EARTH_VELOCITY
        retval[5] /= PyKEP.EARTH_VELOCITY
        retval[6] /= self.__sc.mass
        return retval

    def _compute_constraints_impl(self, x_full):
        import PyKEP
        sc_mass = self.__start_mass
        eqs = []
        ineqs = []

        for i in range(self.__num_legs):
            x = x_full[i * self.__dim_leg:(i + 1) * self.__dim_leg]

            start = PyKEP.epoch(x[0])
            end = PyKEP.epoch(x[0] + x[1])

            # Computing starting spaceraft state
            r, v = self.__seq[i].eph(start)
            x0 = PyKEP.sims_flanagan.sc_state(r, v, sc_mass)

            # Computing ending spaceraft state
            r, v = self.__seq[i + 1].eph(end)
            xe = PyKEP.sims_flanagan.sc_state(r, v, x[3])

            # Building the SF leg
            self.__legs[i].set_spacecraft(PyKEP.sims_flanagan.spacecraft(sc_mass, .3, 3000.))
            self.__legs[i].set(start, x0, x[-3 * self.__nseg:], end, xe)

            # Setting all constraints
            eqs.extend(self.__legs[i].mismatch_constraints())
            ineqs.extend(self.__legs[i].throttles_constraints())

            eqs[-7] /= PyKEP.AU
            eqs[-6] /= PyKEP.AU
            eqs[-5] /= PyKEP.AU
            eqs[-4] /= PyKEP.EARTH_VELOCITY
            eqs[-3] /= PyKEP.EARTH_VELOCITY
            eqs[-2] /= PyKEP.EARTH_VELOCITY
            eqs[-1] /= self.__start_mass

            sc_mass = x[3]  # update mass to final mass of leg

            if i < self.__num_legs - 1:
                x_next = x_full[(i + 1) * self.__dim_leg:(i + 2) * self.__dim_leg]
                time_ineq = x[0] + x[1] + x[2] - x_next[0]
                ineqs.append(time_ineq / 365.25)
            else:
                final_time_ineq = x[0] + x[1] + x[2] - x_full[0] - x_full[-1]  # <- total time
                ineqs.append(final_time_ineq / 365.25)

        retval = eqs + ineqs
        return retval

    def cluster(self, t, with_velocity=True, scaling='astro', eps=0.125, min_samples=10):
        """
        USAGE: cl.cluster(self, t, with_velocity=True, scaling='astro', eps=0.125, min_samples=10)

        - t: epoch (in MJD2000)
        - with_velocity: when True clusters by position and velocity, otherwise only position is used
        - scaling: one of
          - None, or
          - 'standard' (removing mean and scale to standard variance), or
          - 'astro' (scaling by PyKEP.AU and PyKEP.EARTH_VELOCITY)
        - eps: max distance between points in a cluster
        - min_samples: minimum number of samples per cluster
        """
        import PyKEP
        import numpy
        from sklearn.preprocessing import StandardScaler
        from sklearn.cluster import DBSCAN

        self._scaling = scaling
        self._epoch = PyKEP.epoch(t)

        if with_velocity:
            self._X = [
                [elem for tupl in p.eph(self._epoch) for elem in tupl] for p in self._asteroids]
        else:
            self._X = [list(p.eph(self._epoch)[0]) for p in self._asteroids]
        self._X = numpy.array(self._X)

        self._scaler = None
        if self._scaling == 'standard':
            self._scaler = StandardScaler().fit(self._X)
            self._X = self._scaler.transform(self._X)
        elif self._scaling == 'astro':
            scaling_vector = [PyKEP.AU] * 3
            if with_velocity:
                scaling_vector += [PyKEP.EARTH_VELOCITY] * 3
            scaling_vector = numpy.array(scaling_vector)
            self._X = self._X / scaling_vector[None, :]

        self._db = DBSCAN(eps=eps, min_samples=min_samples).fit(self._X)
        self._core_samples = self._db.core_sample_indices_

        self.labels = self._db.labels_
        self.n_clusters = len(
            set(self.labels)) - (1 if -1 in self.labels else 0)

        self.members = {}
        self.core_members = {}
        for label in set(self.labels):
            if int(label) == -1:
                continue
            self.members[int(label)] = [index[0]
                                        for index in numpy.argwhere(self.labels == label)]
            self.core_members[int(label)] = [
                index for index in self._core_samples if self.labels[index] == label]

        if self._scaling == 'standard':
            self._X = self._scaler.inverse_transform(self._X)
        elif self._scaling == 'astro':
            self._X = self._X * scaling_vector[None, :]

def lambert_leg(P1, P2, i, j, t1, t2, tof, vrel=None, dv_launch=0., rendezvous=False):
    """Compute a lambert leg from planet to planet.

    Arguments:
    p1 -- starting planet (str or PyKEP.planet object)
    p2 -- final planet (str or PyKEP.planet object)
    t0 -- start time of leg in MJD2000
    tof -- time of flight in days
    
    Keyword arguments:
    vrel -- caresian coordinates of the relative velocity before the flyby at p1
    dv_launch -- dv discounted at lunch (i.e. if vrel is None)
    rendezvous -- add final dv

    Returns:
    dV, vrel_out, where vrel_out is the relative velocity at the end of the leg at p2
    """

    p1 = PLANETS[str(P1)]
    p2 = PLANETS[str(P2)]

    r1 = state_rosetta.EPH[i][t1][0]
    v1 = state_rosetta.EPH[i][t1][1]
    r2 = state_rosetta.EPH[j][t2][0]
    v2 = state_rosetta.EPH[j][t2][1]

    lambert = kep.lambert_problem(r1, r2, tof * kep.DAY2SEC, p1.mu_central_body, False, 0)

    vrel_in = tuple(map(lambda x, y: x - y, lambert.get_v1()[0], v1))
    vrel_out = tuple(map(lambda x, y: x - y, lambert.get_v2()[0], v2))

    if vrel is None:
        # launch
        dv = max(np.linalg.norm(vrel_in) - dv_launch, 0)
    else:
        # flyby
        #print p1.name, p2.name, np.linalg.norm(vrel_in), np.linalg.norm(vrel_out)
        dv = kep.fb_vel(vrel, vrel_in, p1)

    if rendezvous:
        dv += np.linalg.norm(vrel_out)
        
    return dv, vrel_out

 def pretty(self, x):
     """Decodes the decision vector x"""
     import PyKEP
     start = PyKEP.epoch(x[0])
     end = PyKEP.epoch(x[0] + x[1])
     r, v = self.__departure.eph(start)
     v_list = list(v)
     v_list[0] += x[3]
     v_list[1] += x[4]
     v_list[2] += x[5]
     x0 = PyKEP.sims_flanagan.sc_state(r, v_list, self.__sc.mass)
     r, v = self.__target.eph(end)
     xe = PyKEP.sims_flanagan.sc_state(r, v, x[2])
     self.__leg.set(start, x0, x[-3 * self.__nseg:], end, xe)
     print("A direct interplantary transfer\n")
     print("FROM:")
     print(self.__departure)
     print("TO:")
     print(self.__target)
     print(self.__leg)

    def cluster(self, t, eps=0.125, min_samples=10, metric='orbital', T=180, ref_r=AU, ref_v=EARTH_VELOCITY):
        """
        USAGE: cl.cluster(t, eps=0.125, min_samples=10, metric='orbital', T=180, ref_r=AU, ref_v=EARTH_VELOCITY):

        - t: epoch (in MJD2000)
        - eps: max distance between points in a cluster
        - min_samples: minimum number of samples per cluster
        - metric: one of 'euclidean', 'euclidean_r', orbital'
        - T: average transfer time (used in the definition of the 'orbital' metric)
        - ref_r         reference radius   (used as a scaling factor for r if the metric is 'euclidean' or 'euclidean_r')
        - ref_v         reference velocity (used as a scaling factor for v if the metric is 'euclidean')
        """
        import PyKEP
        import numpy
        from sklearn.cluster import DBSCAN

        self._epoch = PyKEP.epoch(t)

        if metric == 'euclidean':
            self._X = [
                [elem for tupl in p.eph(self._epoch) for elem in tupl] for p in self._asteroids]
            scaling_vector = [ref_r] * 3
            scaling_vector += [ref_v] * 3
        elif metric == 'euclidean_r':
            self._X = [list(p.eph(self._epoch)[0]) for p in self._asteroids]
            scaling_vector = [ref_r] * 3
        elif metric == 'orbital':
            self._T = T
            self._X = [self._orbital_metric(*p.eph(self._epoch)) for p in self._asteroids]
            scaling_vector = [1.] * 6  # no scaling
        self._X = numpy.array(self._X)

        scaling_vector = numpy.array(scaling_vector)
        self._X = self._X / scaling_vector[None, :]

        self._db = DBSCAN(eps=eps, min_samples=min_samples).fit(self._X)
        self._core_samples = self._db.core_sample_indices_

        self.labels = self._db.labels_
        self.n_clusters = len(
            set(self.labels)) - (1 if -1 in self.labels else 0)

        self.members = {}
        self.core_members = {}
        for label in set(self.labels):
            if int(label) == -1:
                continue
            self.members[int(label)] = [index[0]
                                        for index in numpy.argwhere(self.labels == label)]
            self.core_members[int(label)] = [
                index for index in self._core_samples if self.labels[index] == label]

        self._X = self._X * scaling_vector[None, :]

	def __init__(self,mass=1000,Tmax=0.05,Isp=2500,Vinf_0=3,Vinf_f=0,nseg=10,departure = None, target = None, optimise4mass = False):
		"""__init__(self,mass=1000,Tmax=0.05,Isp=2500,Vinf_0=3,Vinf_f=0,nseg=10,departure = erath, target = mars)"""
		try:
			import PyKEP
		except ImportError:
			print("Error while trying 'import PyKEP': is PyKEP installed?")
			raise
		if departure is None:
			departure = PyKEP.planet_ss('earth')
		if target is None:
			target = mars = PyKEP.planet_ss('mars')
		super(py_pl2pl,self).__init__(9 + nseg*3,0,1,9 + nseg,nseg+2,1e-5)
		self.__departure = departure
		self.__target = target
		self.__sc = PyKEP.sims_flanagan.spacecraft(mass,Tmax,Isp)
		self.__Vinf_0 = Vinf_0*1000
		self.__Vinf_f = Vinf_f*1000
		self.__leg = PyKEP.sims_flanagan.leg()
		self.__leg.set_mu(PyKEP.MU_SUN)
		self.__leg.set_spacecraft(self.__sc)
		self.__nseg = nseg
		self.set_bounds([0,10,self.__sc.mass/10,-abs(self.__Vinf_0),-abs(self.__Vinf_0),-abs(self.__Vinf_0),-abs(self.__Vinf_f),-abs(self.__Vinf_f),-abs(self.__Vinf_f)] + [-1] * 3 *nseg,[3000,1500,self.__sc.mass,abs(self.__Vinf_0),abs(self.__Vinf_0),abs(self.__Vinf_0),abs(self.__Vinf_f),abs(self.__Vinf_f),abs(self.__Vinf_f)] + [1] * 3 * nseg)
		self.__optimise4mass = optimise4mass

 def mjd_vts(self, time):
     """Method converting a MJD epoch to the vts format.
     Returns a tuple: the integer part of the MJD and the fractional
     part (converted to seconds) of the MJD. 
     """
     mjd_integer_part = int(pk.epoch(time).mjd // 1)
     mjd_decimal_part = (pk.epoch(time).mjd % 1) * pk.DAY2SEC
     return mjd_integer_part, mjd_decimal_part

    def vts(self, dates):
	# pas compris, si un objet est de type "Date", sa méthode t.vts(...) utilise quoi et retourne quoi?
	"""Method converting an epoch to the vts format.
	it returns the integer part of the MJD and
	the fractional part (converted to seconds) of the MJD. 
	"""
	self.time_integer_part = int(pk.epoch(dates).mjd // 1)
	self.time_decimal_part = (pk.epoch(dates).mjd % 1) * pk.DAY2SEC
	return self.time_integer_part, self.time_decimal_part

def search_min_dv_to_body(body, departure_date, departure_position, host_velocity, min_time_of_flight, time_delta, number):
    time_range = [min_time_of_flight + time_delta*i for i in xrange(number)]
    body_positions = [body.eph(time + departure_date.mjd2000)[0] for time in time_range]
    departure_velocities = [pk.lambert_problem(departure_position, pos, time*pk.DAY2SEC, pk.MU_SUN, False, 0).get_v1()[0] for pos, time in zip(body_positions, time_range)]
    deltaV = [np.linalg.norm(np.array(velocity)-host_velocity) for velocity in departure_velocities]
    index_min = np.array(deltaV).argmin()
    return index_min, departure_velocities[index_min]

def jde_mga_1dsm(seq, t0, tof, slack=5, pop_size=50, n_evolve=10, dv_launch=6127., verbose=False):
    """Runs jDE with mga_1dsm problem."""
    from PyGMO.problem import mga_1dsm_tof
    from PyGMO.algorithm import jde
    from PyGMO import population

    prob = mga_1dsm_tof(seq=[kep.planet_ss(name) for name in seq],
                                   t0=[kep.epoch(t0-slack), kep.epoch(t0+slack)],
                                   tof=[[t-slack, t+slack] for t in tof],
                                   vinf=[0., dv_launch/1000.],
                                   add_vinf_arr=False)
    algo = jde(gen=500, memory=True)
    pop = population(prob, pop_size)
    if verbose:
        print pop.champion.f[0]
    for i in xrange(n_evolve):
        pop = algo.evolve(pop)
        if verbose:
            print pop.champion.f

def lambert_leg(P1, P2, t0, tof):
    
    ast1 = ASTEROIDS[P1]
    ast2 = ASTEROIDS[P2]

    r1, v1 = ast1.eph(kep.epoch(t0))
    r2, v2 = ast2.eph(kep.epoch(t0 + tof))

    lambert = kep.lambert_problem(r1, r2, tof * kep.DAY2SEC, ast1.mu_central_body)

    vrel_in = tuple(map(lambda x, y: -x + y, lambert.get_v1()[0], v1))
    vrel_out = tuple(map(lambda x, y: -x + y, lambert.get_v2()[0], v2))

    dv_lambert = np.linalg.norm(vrel_out) + np.linalg.norm(vrel_in)

    a, _, _, dv_damon = kep.damon(vrel_in, vrel_out, tof*kep.DAY2SEC)
    m_star = kep.max_start_mass(np.linalg.norm(a), dv_damon, T_max, Isp)
    
    return dv_lambert, dv_damon, m_star

def write_output_vts(output_vts_format_file, trajectory, earth, mars):
	""" Function for generating the trajectory of the CubeSat by traj.xyzv 
	in the VTS format """
	def tab_write(value):
        	output_vts_format_file.write('%s\t' % value)

	for value in trajectory:
        	time = value[0]
        	pos = value[1:4]
		time_integer_part = int(pk.epoch(time).mjd // 1) # integer part of mjd time 
		time_decimal_part = (pk.epoch(time).mjd % 1)*pk.DAY2SEC # converting the decimal part of mjd time to seconds
        	tab_write(time_integer_part)
		tab_write(time_decimal_part)
        	tab_write(pos[0]/1000.)  # the position of the CubeSat along the X axis (in km)
        	tab_write(pos[1]/1000.)
        	tab_write(pos[2]/1000.)
        	tab_write(value[4])      # the velocity of the CubeSat along the X axis (in m/s)
        	tab_write(value[5])
       		tab_write(value[6])
        	output_vts_format_file.write('\n')