From 4ea9d5149d3f561f1a29689d716d7b446a742e64 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Wed, 16 Jan 2013 19:12:50 -0600 Subject: [PATCH 01/25] Add a parameter for the shape of the gravity turn. --- ascent.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/ascent.py b/ascent.py index 23da85f..ab4b792 100644 --- a/ascent.py +++ b/ascent.py @@ -63,7 +63,8 @@ def __init__(self, initialVelocity = None, launchInclination = 0, acceleration = None, - timestep = 1): + timestep = 1, + gravityTurnCurve = 1): """ Compute a climb slope for exiting the atmosphere and achieving orbit. @@ -199,6 +200,7 @@ def __str__(self): else: # What fraction of the turn have we done? ratio = (alt - gravityTurnStart) / (gravityTurnEnd - gravityTurnStart) + ratio **= gravityTurnCurve # The more we did, the closer we want to be to pi/2 from # the straight-up orientation. phiSurf = - (ratio * math.pi / 2) From df290602d9390def8999f6dbd32ffa0add198279 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Wed, 16 Jan 2013 19:16:44 -0600 Subject: [PATCH 02/25] Add apoapsis and periapsis values. --- planet.py | 46 ++++++++++++++++++++++++++-------------------- 1 file changed, 26 insertions(+), 20 deletions(-) diff --git a/planet.py b/planet.py index ba486ae..36b66a1 100644 --- a/planet.py +++ b/planet.py @@ -39,7 +39,7 @@ class planet(object): - def __init__(self, name, gravityParam, SOI, radiusKm, siderealPeriod, datumPressure, scale): + def __init__(self, name, gravityParam, SOI, radiusKm, siderealPeriod, datumPressure, scale, ap, pe): self.name = name self.mu = gravityParam # m^3/s^2 self.SOI = SOI # m @@ -49,6 +49,9 @@ def __init__(self, name, gravityParam, SOI, radiusKm, siderealPeriod, datumPress self.siderealRotationSpeed = 2 * pi * self.radius / siderealPeriod # m/s self.datumPressure = datumPressure # atm self.scale = scale # m (pressure falls by e every scale altitude) + self.ap = ap + self.pe = pe + self.sma = (ap + pe) / 2 def __str__(self): return self.name @@ -295,27 +298,30 @@ def topOfAtmosphere(self): # the sun has an infinite SOI; arbitrarily set it to 500x the orbit of Jool sunSOI = 500 * 71950638386 +# the sun doesn't orbit anything +sunAp = sunPe = 0 + planets = dict([ (p.name.lower(), p) for p in ( - # name, mu, SOI, radiusKm, sidereal, atm, scale - planet("Kerbol", 1.172332794E+18, sunSOI, 261600, 432000, 0, 0), - planet("Moho", 245250003655, 11206449 , 250, 1210000, 0, 0), - planet("Eve", 8171730229211, 85109365 , 700, 80500, 5, 7000), - planet("Gilly", 8289450, 126123.27, 13, 28255, 0, 0), - planet("Kerbin", 3.5316E+12, 84159286 , 600, 21600, 1, 5000), - planet("Mun", 65138397521, 2429559.1 , 200, 138984, 0, 0), - planet("Minmus", 1765800026, 2247428.4 , 60, 40400, 0, 0), - planet("Duna", 301363211975, 47921949 , 320, 65518, 0.2, 3000), - planet("Ike", 18568368573, 1049598.9 , 130, 65518, 0, 0), - planet("Dres", 21484488600, 32700000 , 138, 34800, 0, 0), - planet("Jool", 282528004209995, 2455985200 , 600, 36000, 15, 9000), - planet("Laythe", 1962000029236, 3723645.8 , 500, 52981, 0.8, 4000), - planet("Vall", 207481499474, 2406401.4 , 300, 105962, 0, 0), - planet("Tylo", 2825280042100, 10856518 , 600, 211926, 0, 0), - planet("Bop", 2486834944, 993002.8 , 65, 544507, 0, 0), - planet("Pol", 227905920, 1041613 , 44, 901903, 0, 0), - planet("Eeloo", 74410814527, 119087000 , 210, 19460, 0, 0), + # name, mu, SOI, radiusKm, sidereal, atm, scale, apoapsis, periapsis + planet("Kerbol", 1.172332794E+18, sunSOI, 261600, 432000, 0, 0, sunAp, sunPe), + planet("Moho", 245250003655, 11206449 , 250, 1210000, 0, 0, 6315765980, 4210510628), + planet("Eve", 8171730229211, 85109365 , 700, 80500, 5, 7000, 9931011387, 9734357701), + planet("Gilly", 8289450, 126123.27, 13, 28255, 0, 0, 48825000, 14175000), + planet("Kerbin", 3.5316E+12, 84159286 , 600, 21600, 1, 5000, 13599840256, 13599840256), + planet("Mun", 65138397521, 2429559.1 , 200, 138984, 0, 0, 12000000, 12000000), + planet("Minmus", 1765800026, 2247428.4 , 60, 40400, 0, 0, 47000000, 47000000), + planet("Duna", 301363211975, 47921949 , 320, 65518, 0.2, 3000, 21783189163, 19669121365), + planet("Ike", 18568368573, 1049598.9 , 130, 65518, 0, 0, 3296000, 3104000), + planet("Dres", 21484488600, 32700000 , 138, 34800, 0, 0, 46761053522, 34917642884), + planet("Jool", 282528004209995, 2455985200 , 600, 36000, 15, 9000, 71950638386, 65073282253), + planet("Laythe", 1962000029236, 3723645.8 , 500, 52981, 0.8, 4000, 27184000, 27184000), + planet("Vall", 207481499474, 2406401.4 , 300, 105962, 0, 0, 43152000, 43152000), + planet("Tylo", 2825280042100, 10856518 , 600, 211926, 0, 0, 68500000, 68500000), + planet("Bop", 2486834944, 993002.8 , 65, 544507, 0, 0, 129057500, 79942500), + planet("Pol", 720792000, 1041613 , 44, 901903, 0, 0, 210624206, 149155794), + planet("Eeloo", 74410814527, 119087000 , 210, 19460, 0, 0, 113549713200, 66687926800), ) ]) -del sunSOI +del sunSOI, sunAp, sunPe def getPlanet(name): return planets[name.lower()] From c94a76893d202da1dbbecb55bfeb8eb92dc39180 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Wed, 16 Jan 2013 20:14:32 -0600 Subject: [PATCH 03/25] Allow the coefficient of drag to be specified. --- ascent.py | 10 +++++++--- planet.py | 37 +++++++++++++++++++------------------ 2 files changed, 26 insertions(+), 21 deletions(-) diff --git a/ascent.py b/ascent.py index ab4b792..04e537a 100644 --- a/ascent.py +++ b/ascent.py @@ -64,7 +64,8 @@ def __init__(self, launchInclination = 0, acceleration = None, timestep = 1, - gravityTurnCurve = 1): + gravityTurnCurve = 1, + dragCoefficient = None): """ Compute a climb slope for exiting the atmosphere and achieving orbit. @@ -141,6 +142,9 @@ def __str__(self): else: v = list(initialVelocity) + if dragCoefficient is None: + dragCoefficient = planet.defaultDragCoefficient + # p and v are the current position and velocity; we update these. # We use a coordinate system centered at the core of the planet, # aligned so that the initial altitude forms the y axis. @@ -209,7 +213,7 @@ def __str__(self): # Compute the gravity and drag. a_grav = planet.gravity(alt) - a_drag = planet.drag(alt, L2(v)) + a_drag = planet.drag(alt, L2(v), dragCoefficient) # Terminal velocity for our altitude: #print ("drag %g m/s^2, grav %g m/s^2, term %g m/s" % (a_drag, a_grav, v_term)) @@ -250,7 +254,7 @@ def thrustResult(thrust): def findTerminalThrust(): # above the top of the atmosphere, there is no limit! if alt >= TOA: return 1e30 - v_term = planet.terminalVelocity(alt) + v_term = planet.terminalVelocity(alt, dragCoefficient) a = timestep * timestep v_noTX = v_nothrust[0] * cos(phi) + v_nothrust[1] * sin(phi) b = 2 * timestep * v_noTX diff --git a/planet.py b/planet.py index 36b66a1..69969db 100644 --- a/planet.py +++ b/planet.py @@ -25,16 +25,13 @@ """ # The drag model in KSP 0.18.1 is -# F_{drag} = gamma D v^2 m -# where D is atmospheric density, v is the speed, m the mass (as a standin -# for cross-section). gamma combines a bunch of coefficients in one -# (including the 1/2 that would normally be there). It is not actually -# a constant, but it varies little over reported terminal velocities and over -# the parts we actually use (though aero components have lower drag). -coefficientOfDrag = 0.2 # Assumed constant, not really the case +# F_{drag} = gamma D P v^2 m +# where D is the coefficient of drag, P is atmospheric pressure, v is the speed, +# m the mass (as a standin for cross-section). gamma combines a bunch of +# coefficients in one (including the 1/2 that would normally be there). dragMultiplier = 0.008 kerbinSurfaceDensity = 1.2230948554874 -gamma = 0.5 * coefficientOfDrag * dragMultiplier * kerbinSurfaceDensity +gamma = 0.5 * dragMultiplier * kerbinSurfaceDensity class planet(object): @@ -52,6 +49,7 @@ def __init__(self, name, gravityParam, SOI, radiusKm, siderealPeriod, datumPress self.ap = ap self.pe = pe self.sma = (ap + pe) / 2 + self.defaultDragCoefficient = 0.2 def __str__(self): return self.name @@ -229,7 +227,7 @@ def soiBurn(self, altitude, soiSpeed): return v - self.orbitalVelocity(altitude) - def terminalVelocity(self, altitude): + def terminalVelocity(self, altitude, dragCoefficient = None): """ Return the terminal velocity at a given altitude. This is the speed at which drag is as strong as gravity. @@ -237,14 +235,15 @@ def terminalVelocity(self, altitude): altitude is in meters. """ - # v = sqrt(g/ (datumPressure * gamma * e^{-altitude/scale})) + # v = sqrt(g/ (datumPressure * gamma * D * e^{-altitude/scale})) # Then pull the e term out, and remember that gravity changes # (slightly) with altitude. + if dragCoefficient is None: dragCoefficient = self.defaultDragCoefficient if altitude is None or altitude >= self.topOfAtmosphere(): return float("inf") return exp(0.5 * altitude / self.scale) * sqrt( - self.gravity(altitude) / (gamma * self.datumPressure) ) + self.gravity(altitude) / (gamma * dragCoefficient * self.datumPressure) ) - def drag(self, altitude, velocity): + def drag(self, altitude, velocity, dragCoefficient = None): """ Return the acceleration due to drag, in m/s^2. @@ -253,13 +252,15 @@ def drag(self, altitude, velocity): altitude is in meters, velocity in m/s """ - # F_{drag} = D v^2 m gamma, where D is atmospheric pressure, v is - # the speed, m the mass (as a standin for cross-section), and gamma - # combines a bunch of variables that are close enough to constant. - # To get the acceleration, divide by mass, which cancels out m: - # a_{drag} = gamma P v^2 + # F_{drag} = P v^2 m D gamma, where P is atmospheric pressure, v is + # the speed, m the mass (as a standin for cross-section), D is the + # coefficient of drag, and gamma combines a bunch of variables that are + # close enough to constant. To get the acceleration, divide by mass, + # which cancels out m: + # a_{drag} = gamma P v^2 D + if dragCoefficient is None: dragCoefficient = self.defaultDragCoefficient if altitude is None or altitude >= self.topOfAtmosphere(): return 0 - return gamma * self.pressure(altitude) * velocity * velocity + return gamma * dragCoefficient * self.pressure(altitude) * (velocity ** 2) def pressure(self, altitude): """ From 841d69db48da12605c8ba5eb602e5739a01abc58 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Fri, 1 Feb 2013 17:56:40 -0600 Subject: [PATCH 04/25] Add optional apo/peri args to orbitalVelocity(). --- planet.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/planet.py b/planet.py index 69969db..a567a63 100644 --- a/planet.py +++ b/planet.py @@ -63,7 +63,7 @@ def gravity(self, altitude = 0): r = self.radius + altitude return self.mu / (r*r) - def orbitalVelocity(self, altitude): + def orbitalVelocity(self, altitude, ap = None, pe = None): """ Return the velocity required to maintain an orbit with the given altitude at the semimajor axis (e.g. a circular orbit at that @@ -71,7 +71,11 @@ def orbitalVelocity(self, altitude): altitude is in meters. """ - return sqrt(self.mu / (self.radius + altitude)) + if ap is None: ap = altitude + if pe is None: pe = altitude + r = self.radius + altitude + a = self.radius + (ap + pe) / 2 + return sqrt(self.mu * (2 / r - 1 / a)) def escapeVelocity(self, altitude): """ From d78a875a693087e4724088068b4ddfaf13d6ba00 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Wed, 20 Feb 2013 21:28:20 -0600 Subject: [PATCH 05/25] Learn optimal ascents with a genetic algorithm. --- learn-ascent.py | 227 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 227 insertions(+) create mode 100644 learn-ascent.py diff --git a/learn-ascent.py b/learn-ascent.py new file mode 100644 index 0000000..f60cdda --- /dev/null +++ b/learn-ascent.py @@ -0,0 +1,227 @@ +import argparse +import math +import os +import random +import sys + +import ascent +import planet + +class Profile: + + # Class variables. + planet = None + alt0 = None + alt1 = None + accel = None + drag = None + + def __init__(self, gt0, gt1, curve): + self.gt0 = sorted([0, gt0, self.alt1])[1] + self.gt1 = sorted([0, gt1, self.alt1])[1] + self.curve = sorted([0, curve, 2])[1] + self._calc_score() + self.generation = 1 + + @classmethod + def init(cls, planet, alt0, alt1, accel, drag): + cls.planet = planet + cls.alt0 = alt0 + cls.alt1 = alt1 + cls.accel = accel + cls.drag = drag + + @classmethod + def from_string(cls, str): + tokens = [float(f) for f in str.strip().split(" ")] + if len(tokens) == 4: + tokens.append(-1) + (gt0, gt1, curve, score) = tokens[:4] + return Profile(gt0, gt1, curve) + + #9.72 23.87 0.3381 1119.067026 + @classmethod + def random(cls): + gt0 = random.random() * (cls.alt1 / 2) + gt1 = gt0 + random.random() * (cls.alt1 - gt0) + curve = random.random() * 2 + return Profile(gt0, gt1, curve) + + def mutated(self): + return Profile( self._mutate(self.gt0), + self._mutate(self.gt1), + self._mutate(self.curve)) + + def _mutate(self, value): + amount = 0.01 + + # Mutate more slowly in later generations. + #amount /= math.log(self.generation + 1) + + return value + (random.random() - 0.5) * math.sqrt(value) * amount + + def _calc_score(self): + try: + self.ascent = ascent.climbSlope(self.planet, + orbitAltitude = self.alt1 * 1000, + gravityTurnStart = self.gt0 * 1000, + gravityTurnEnd = self.gt1 * 1000, + gravityTurnCurve = self.curve, + acceleration = self.planet.gravity() * self.accel, + initialAltitude = self.alt0 * 1000, + dragCoefficient = self.drag + ) + self.score = self.ascent.deltaV() + except ascent.BadFlightPlanException as bfpe: + self.score = -1 + self.ascent = None + + def _combine(self, a, b): + average = (a ** 2 + b ** 2) ** 0.5 + return self._mutate(average) + + def __lt__(self, other): + return self.score != -1 and self.score < other.score + + def combine(self, other): + return Profile(self._combine(self.gt0, other.gt0), + self._combine(self.gt1, other.gt1), + self._combine(self.curve, other.curve)) + + def better_than(self, other): + return (not other) or other.score < 0 or (self.score > 0 and self.score < other.score) + + def worse_than(self, other): + return (not other) or self.score > other.score + + def __str__(self): + return "%f %f %f %f" % (self.gt0, self.gt1, self.curve, self.score) + + def guide(self): + angleStep = 15 + guide = "(angle alt) " + for angle in range(0, 90 + 1, angleStep): + dy = self.gt1 - self.gt0 + alt = self.gt0 + ((float(angle) / 90) ** (1.0 / self.curve)) * dy + if angle != 0: + guide += ", " + guide += "%d %4.1f" % (angle, alt) + guide += "\n" + if self.ascent: + guide += "average atmospheric pressure: %f\n" % (self.ascent.avgAtm) + guide += "dv to circularize at apoapsis: %f\n" % (self.ascent.dV_circ) + #guide += str(self.ascent) + "\n" + return guide + +def select(pool, s): + for profile in pool: + if random.random() > s: + return profile + return pool[-1] + + +def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, poolSize = 20, fileIn = None): + + p = planet.planets[planetName.lower()] + + if endAlt is None: + endAlt = math.ceil(p.topOfAtmosphere() / 5000) * 5 + + print("ascending to %d km with %.2f g's of acceleration" % (endAlt, accel)) + Profile.init(p, startAlt, endAlt, accel, drag) + + fileOut = "%s_%d_%d_%.2f_%.2f.txt" % (p.name, startAlt, endAlt, accel, drag) + + if fileIn is None: + fileIn = fileOut + + pool = [] + + if fileIn and os.path.exists(fileIn): + with open(fileIn) as f: + for line in f.readlines(): + pool.append(Profile.from_string(line)) + + while len(pool) < poolSize: + profile = Profile.random() + pool.append(profile) + + bestEver = None + gen = 1 + needNewline = False + + try: + while True: + best = None + worst = None + total = 0 + successes = 0 + # 7.6 45.0 0.653 = 4435.43 + # 8.7 45.9 0.610 = 4436.00 + candidates = [] + for profile in pool: + profile.generation = gen + if profile.better_than(best): + best = profile + if profile.better_than(bestEver): + bestEver = profile + if (needNewline): + sys.stdout.write("\n") + needNewline = False + print(" " * 8 + str(profile)) + print(profile.guide()) + + if profile.worse_than(worst): + worst = profile + + if profile.score > 0: + total += profile.score + successes += 1 + candidates.append(profile) + candidates.sort() + if successes: + needNewline = True + sys.stdout.write("\r%6d: best %f, average %f" % (gen, best.score, total / successes)) + newPool = [] + SELECT_P = 0.5 + + # Automatically select the top candidate and a mutant of it. + if candidates: + newPool.append(candidates[0]) + newPool.append(candidates[0].mutated()) + + while len(newPool) < min(poolSize / 2, successes): + a = select(candidates, SELECT_P) + b = select(candidates, SELECT_P) + newPool.append(a.combine(b)) + + while len(newPool) < poolSize: + newPool.append(Profile.random()) + pool = newPool + gen += 1 + + except KeyboardInterrupt: + pool.sort() + with open(fileOut, "w") as f: + for profile in pool: + f.write("%s\n" % profile) + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description = "Learn an ascent.") + parser.add_argument('planet', metavar='planet', + help='planet or moon') + parser.add_argument('alt0', metavar='alt0', type=int, + help='initial altitude (km; defaults to 0)') + parser.add_argument('alt1', metavar='alt1', type=int, + help='final altitude (km; defaults to the nearest 5 km mark above the atmosphere)') + parser.add_argument('-a', metavar='accel', type=float, default=2, + help='ship acceleration as a multiple of planet surface gravity (default: %(default)s)') + parser.add_argument('-d', metavar='drag', type=float, default=0.2, + help='drag coefficient (default: %(default)s)') + parser.add_argument('-p', metavar='poolSize', type=int, default=20, + help='pool size (default: %(default)s)') + #parser.print_help() + + args = parser.parse_args(sys.argv[1:]) + + learnAscent(args.planet, args.alt0, args.alt1, args.a, args.d, args.p) From d128e63b421f8c5299f4d93fcf6ce9f37af7ed17 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Wed, 20 Feb 2013 23:30:46 -0600 Subject: [PATCH 06/25] Allow variable acceleration to be modeled. --- ascent.py | 48 +++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 39 insertions(+), 9 deletions(-) diff --git a/ascent.py b/ascent.py index 04e537a..dc3b5d9 100644 --- a/ascent.py +++ b/ascent.py @@ -65,7 +65,9 @@ def __init__(self, acceleration = None, timestep = 1, gravityTurnCurve = 1, - dragCoefficient = None): + dragCoefficient = None, + specificImpulse = None, + shipThrust = None): """ Compute a climb slope for exiting the atmosphere and achieving orbit. @@ -83,11 +85,14 @@ def __init__(self, with an angle that varies linearly with altitude. The start and end points of the turn can be changed, the curve cannot. - There's no model of variable acceleration. It would be easy - to add in since we are using numerical integration; you'd want to - take in a function from (deltaV, altitude, velocity) -> - acceleration to model jets. By default, acceleration is assumed to be 2.2g, - which is good low down (except the first ~100m) but too little at altitude. + Variable acceleration is modeled if shipThrust (in kN) and specificImpulse + (in s) are specified. Initial ship mass is calculated by assuming that the given + acceleration is the maximum acceleration achievable at launch. At each step of + the ascent, the total delta-V achieved and the given specificImpulse are used to + reduce the ship's mass and thus increase its maximum acceleration. + + By default, acceleration is assumed to be 2.2g, which is good low down (except + the first ~100m) but too little at altitude. Specify the timestep to change the accuracy; by default we use a timestep of 1s. It is a fixed timestep. Smaller timesteps @@ -155,9 +160,13 @@ def __str__(self): t = 0 # s dragLoss = 0 # m/s if acceleration is None: - a_thrust = 2.2 * planet.gravity(initialAltitude) # m/s^2 - else: - a_thrust = acceleration + acceleration = 2.2 * planet.gravity(initialAltitude) # m/s^2 + + variable_accel = not (specificImpulse is None and shipThrust is None) + + if variable_accel: + # f=ma, so m=f/a + mass = shipThrust / acceleration TOA = planet.topOfAtmosphere() # m if orbitAltitude is None: @@ -306,6 +315,27 @@ def findApoapsisThrust(thrust_term): thrust_apo = findApoapsisThrust(thrust_term) thrust = min(thrust_term, thrust_apo) + + def getAcceleration(dv, alt): + acc = acceleration + if variable_accel: + # dv = isp * g0 * ln(m0/m1) + # m0 / m1 = math.exp(dv / (isp * g0)) + # 1 + dm / m0 = math.exp(dv / (isp * g0)) + # dm = m0 * (math.exp(dv / (isp * g0)) - 1) + # m = m0 - dm + # m = m0 * (1 - (math.exp(dv / (isp * g0)) - 1)) + # m = m0 * (2 - math.exp(dv / (isp * g0))) + # m0 = shipThrust / acceleration + # m = (2 - math.exp(dv / (isp * g0))) * shipThrust / acceleration + # acc = shipThrust / m + # acc = acceleration / (2 - math.exp(dv / (isp * g0))) + acc /= (2 - math.exp(dv / (specificImpulse * 9.81))) + return acc + + # f=ma, so a=f/m + a_thrust = getAcceleration(dV, alt) + #print ("thrust: %g to reach term, %g to reach apo" % (thrust_term, thrust_apo)) if thrust < a_thrust: thrustLimit = 0 From a2901c4444b37e44f194d4396bde6da4f1256b6d Mon Sep 17 00:00:00 2001 From: gmcnew Date: Sat, 23 Feb 2013 22:47:49 -0600 Subject: [PATCH 07/25] Make ascent simulations about 40% faster. --- ascent.py | 36 +++++++++++++++++++++++++++--------- planet.py | 23 +++++++++++++---------- 2 files changed, 40 insertions(+), 19 deletions(-) diff --git a/ascent.py b/ascent.py index dc3b5d9..9e5c2f7 100644 --- a/ascent.py +++ b/ascent.py @@ -190,6 +190,11 @@ def __str__(self): # decide how much to burn and update our position. targetApoapsis = orbitAltitude + planet.radius climbSlope = [] # ClimbPoint list + + # Keep track of approximately how much thrust would be needed to reach + # the desired apoapsis. + thrust_apo = None + while h < targetApoapsis and t < 1000: if h < planet.radius: # We crashed! Bring more thrust next time. @@ -234,11 +239,14 @@ def __str__(self): - (a_drag * sin(theta) + a_grav * sin(psi)) ) v_nothrust = [ v[i] + loss[i] * timestep for i in range(2) ] + + def thrustResult2(thrust, tx, ty): + return (v_nothrust[0] + thrust * tx, v_nothrust[1] + thrust * ty) + def thrustResult(thrust): - return ( - v_nothrust[0] + thrust * cos(phi) * timestep, - v_nothrust[1] + thrust * sin(phi) * timestep, - ) + tx = cos(phi) * timestep + ty = sin(phi) * timestep + return thrustResult2(thrust, tx, ty) # Compute the acceleration that gets us to terminal velocity in the # direction of thrust. @@ -274,7 +282,7 @@ def findTerminalThrust(): # Binary search the thrust to achieve the apoapsis. # I'm sure this could be worked out analytically. # Use the terminal velocity thrust to reduce our search space. - def findApoapsisThrust(thrust_term): + def findApoapsisThrust(thrust_term, guess = None): # The most we could want to speed up is enough to immediately # get our momentum to match what it will when we have a @@ -295,14 +303,23 @@ def findApoapsisThrust(thrust_term): (p[0]*cos(thetaSurf) + p[1]*sin(thetaSurf))) thrust_targetMax = (v_targetMax - L2(v_nothrust)) / timestep + tx = cos(phi) * timestep + ty = sin(phi) * timestep + + theta = math.atan2(p[1], p[0]) + thrust_lo = 0 thrust_hi = min(thrust_targetMax, thrust_term) # while we are off by more than 1mm/s, binary search # we waste a bit of time searching for really big solutions. while thrust_hi - thrust_lo > 1e-3: - thrust = (thrust_hi + thrust_lo) / 2 - v_next = thrustResult(thrust) - (apoapsis, _) = planet.determineOrbit2(p, v_next) + if guess and guess < thrust_hi: + thrust = guess + guess = None + else: + thrust = (thrust_hi + thrust_lo) / 2 + v_next = thrustResult2(thrust, tx, ty) + (apoapsis, _) = planet.determineOrbit3(h, theta, v_next) # if apoapsis is too high, or negative (i.e. hyperbolic), reduce thrust if apoapsis > targetApoapsis or apoapsis < 0: @@ -312,7 +329,8 @@ def findApoapsisThrust(thrust_term): return thrust_lo # or hi, it's only 1mm/s difference thrust_term = findTerminalThrust() - thrust_apo = findApoapsisThrust(thrust_term) + guess = thrust_apo + thrust_apo = findApoapsisThrust(thrust_term, guess) thrust = min(thrust_term, thrust_apo) diff --git a/planet.py b/planet.py index a567a63..457e5c3 100644 --- a/planet.py +++ b/planet.py @@ -126,28 +126,31 @@ def determineOrbit(self, h, velocity): p = h * v * cos(theta) return quadratic(ainv, -2, p*p / self.mu) - def determineOrbit2(self, position, velocity): - """ - Given a position and velocity in some coordinate frame centered about - the center of the planet (doesn't matter, as long as 1 unit is 1 - meter), return the apsides, in meters above the core. - If one of them is negative, the orbit is hyperbolic. - """ - + def determineOrbit3(self, h, theta, velocity): + # h is distance from the planet's core. # theta is the angle of the position above the horizontal. # phi is the angle of the velocity above the horizontal. # psi is the angle of the velocity above the surface tangent. # psi = (90 - theta) + phi # thetabar = 90-theta - theta = math.atan2(position[1], position[0]) phi = math.atan2(velocity[1], velocity[0]) thetabar = (math.pi/2) - theta psi = thetabar + phi - h = L2(position) v = L2(velocity) return self.determineOrbit(h, (v, psi)) + def determineOrbit2(self, position, velocity): + """ + Given a position and velocity in some coordinate frame centered about + the center of the planet (doesn't matter, as long as 1 unit is 1 + meter), return the apsides, in meters above the core. + If one of them is negative, the orbit is hyperbolic. + """ + theta = math.atan2(position[1], position[0]) + h = L2(position) + return self.determineOrbit3(h, theta, velocity) + def hohmann(self, a1, a2): """ Given a circular orbit at altitude a1 and a target orbit at altitude From f3120ee3965ff525d40be0b49afea3c507475443 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Sun, 24 Feb 2013 15:16:35 -0600 Subject: [PATCH 08/25] Make L2() about 15% faster by replacing sum(). --- physics.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/physics.py b/physics.py index abccefa..20c4464 100644 --- a/physics.py +++ b/physics.py @@ -34,4 +34,7 @@ def quadratic(a, b, c): return ( (-b + sqrtdiscriminant) / a2, (-b - sqrtdiscriminant) / a2) def L2(vector): - return math.sqrt(sum(x*x for x in vector)) + a = 0 + for x in vector: + a += x ** 2 + return math.sqrt(a) From 77cbee4b0a612a58cd4a424dbd450bbc9c7449ff Mon Sep 17 00:00:00 2001 From: gmcnew Date: Sun, 24 Feb 2013 15:41:13 -0600 Subject: [PATCH 09/25] Increase speed 60%: limit thrust_apo to a_thrust. --- ascent.py | 34 ++++++++++++++++++++-------------- 1 file changed, 20 insertions(+), 14 deletions(-) diff --git a/ascent.py b/ascent.py index 9e5c2f7..921e158 100644 --- a/ascent.py +++ b/ascent.py @@ -65,6 +65,7 @@ def __init__(self, acceleration = None, timestep = 1, gravityTurnCurve = 1, + calculateExtraThrust = False, dragCoefficient = None, specificImpulse = None, shipThrust = None): @@ -135,8 +136,10 @@ def __str__(self): self.deltaV, self.dragLoss, self.thrust, - "" if self.thrustLimited == 0 else - (" needs %g m/s^2 more thrust" % self.thrustLimited) + "" if not self.thrustLimited else + (" needs %smore thrust" % ( + "" if self.thrustLimited is True else + ("%g m/s^2" % self.thrustLimited))) )) self.planet = planet @@ -164,10 +167,6 @@ def __str__(self): variable_accel = not (specificImpulse is None and shipThrust is None) - if variable_accel: - # f=ma, so m=f/a - mass = shipThrust / acceleration - TOA = planet.topOfAtmosphere() # m if orbitAltitude is None: orbitAltitude = TOA + 1000 @@ -328,12 +327,6 @@ def findApoapsisThrust(thrust_term, guess = None): thrust_lo = thrust return thrust_lo # or hi, it's only 1mm/s difference - thrust_term = findTerminalThrust() - guess = thrust_apo - thrust_apo = findApoapsisThrust(thrust_term, guess) - - thrust = min(thrust_term, thrust_apo) - def getAcceleration(dv, alt): acc = acceleration if variable_accel: @@ -351,16 +344,29 @@ def getAcceleration(dv, alt): acc /= (2 - math.exp(dv / (specificImpulse * 9.81))) return acc - # f=ma, so a=f/m a_thrust = getAcceleration(dV, alt) + thrust_term = findTerminalThrust() + guess = thrust_apo + thrust_limit = thrust_term + if not calculateExtraThrust: + # It's nice to know if thrust_apo exceeds possible thrust + # (a_thrust). However, determining thrust_apo precisely is very + # time-consuming, so let's limit it to 1 m/s^2 greater than + # a_thrust. This will still tell us if we need more thrust, it + # just won't tell us how much extra thrust we need. + thrust_limit = min(thrust_term, a_thrust + 1) + thrust_apo = findApoapsisThrust(thrust_limit, guess) + + thrust = min(thrust_term, thrust_apo) + #print ("thrust: %g to reach term, %g to reach apo" % (thrust_term, thrust_apo)) if thrust < a_thrust: thrustLimit = 0 if thrust < 0: thrust = 0 else: - thrustLimit = thrust - a_thrust + thrustLimit = (thrust - a_thrust) if calculateExtraThrust else True thrust = a_thrust # Update everything! From d9f5374e5c07dd43297b914c4bbe15bfe69e963e Mon Sep 17 00:00:00 2001 From: gmcnew Date: Thu, 28 Feb 2013 20:42:43 -0600 Subject: [PATCH 10/25] Miscellaneous tweaks. (Cleanup is likely needed.) --- learn-ascent.py | 83 ++++++++++++++++++++++++++++++++++++------------- 1 file changed, 62 insertions(+), 21 deletions(-) diff --git a/learn-ascent.py b/learn-ascent.py index f60cdda..7727ea7 100644 --- a/learn-ascent.py +++ b/learn-ascent.py @@ -16,11 +16,24 @@ class Profile: accel = None drag = None + cache = {} + MAX_CACHE_SIZE = 10000 + def __init__(self, gt0, gt1, curve): - self.gt0 = sorted([0, gt0, self.alt1])[1] - self.gt1 = sorted([0, gt1, self.alt1])[1] - self.curve = sorted([0, curve, 2])[1] - self._calc_score() + # Limit values to 3 decimal places of precision. + (gt0, gt1, curve) = [(int(x * 1000) / 1000.0) for x in (gt0, gt1, curve)] + + self.gt0 = sorted([0, gt0, self.alt1])[1] + self.gt1 = sorted([0, gt1, self.alt1])[1] + self.curve = sorted([0, curve, 1])[1] + + self.score = self.cache.get((self.gt0, self.gt1, self.curve), None) + if self.score is None: + if len(self.cache) == self.MAX_CACHE_SIZE: + del self.cache[random.choice(self.cache.keys())] + self._calc_score() + self.cache[(self.gt0, self.gt1, self.curve)] = self.score + self.generation = 1 @classmethod @@ -31,6 +44,10 @@ def init(cls, planet, alt0, alt1, accel, drag): cls.accel = accel cls.drag = drag + @classmethod + def clear_cache(cls): + cls.cache = {} + @classmethod def from_string(cls, str): tokens = [float(f) for f in str.strip().split(" ")] @@ -48,12 +65,13 @@ def random(cls): return Profile(gt0, gt1, curve) def mutated(self): - return Profile( self._mutate(self.gt0), - self._mutate(self.gt1), - self._mutate(self.curve)) + vals = [self.gt0, self.gt1, self.curve] + i = random.randint(0, 2) + vals[i] = self._mutate(vals[i]) + return Profile(vals[0], vals[1], vals[2]) def _mutate(self, value): - amount = 0.01 + amount = 0.1 # Mutate more slowly in later generations. #amount /= math.log(self.generation + 1) @@ -95,7 +113,7 @@ def worse_than(self, other): return (not other) or self.score > other.score def __str__(self): - return "%f %f %f %f" % (self.gt0, self.gt1, self.curve, self.score) + return "%.3f %.3f %.3f %f" % (self.gt0, self.gt1, self.curve, self.score) def guide(self): angleStep = 15 @@ -107,10 +125,6 @@ def guide(self): guide += ", " guide += "%d %4.1f" % (angle, alt) guide += "\n" - if self.ascent: - guide += "average atmospheric pressure: %f\n" % (self.ascent.avgAtm) - guide += "dv to circularize at apoapsis: %f\n" % (self.ascent.dV_circ) - #guide += str(self.ascent) + "\n" return guide def select(pool, s): @@ -119,15 +133,17 @@ def select(pool, s): return profile return pool[-1] +SILENT = True -def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, poolSize = 20, fileIn = None): +def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, poolSize = 20, fileIn = None, genLimit = None): p = planet.planets[planetName.lower()] if endAlt is None: endAlt = math.ceil(p.topOfAtmosphere() / 5000) * 5 - print("ascending to %d km with %.2f g's of acceleration" % (endAlt, accel)) + if not SILENT: + print("ascending to %d km with %.2f g's of acceleration" % (endAlt, accel)) Profile.init(p, startAlt, endAlt, accel, drag) fileOut = "%s_%d_%d_%.2f_%.2f.txt" % (p.name, startAlt, endAlt, accel, drag) @@ -137,7 +153,7 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, pool = [] - if fileIn and os.path.exists(fileIn): + if fileIn and os.path.exists(fileIn) and not SILENT: with open(fileIn) as f: for line in f.readlines(): pool.append(Profile.from_string(line)) @@ -148,8 +164,10 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, bestEver = None gen = 1 + lastChange = 0 needNewline = False + bestThisRound = None try: while True: best = None @@ -161,10 +179,13 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, candidates = [] for profile in pool: profile.generation = gen - if profile.better_than(best): - best = profile + if profile.better_than(best) and not SILENT: if profile.better_than(bestEver): bestEver = profile + best = profile + if profile.better_than(bestThisRound): + lastChange = gen + bestThisRound = profile if (needNewline): sys.stdout.write("\n") needNewline = False @@ -181,7 +202,8 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, candidates.sort() if successes: needNewline = True - sys.stdout.write("\r%6d: best %f, average %f" % (gen, best.score, total / successes)) + if not SILENT: + sys.stdout.write("\r%6d: best %f, average %f, best ever %f" % (gen, best.score, total / successes, bestEver.score)) newPool = [] SELECT_P = 0.5 @@ -190,6 +212,13 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, newPool.append(candidates[0]) newPool.append(candidates[0].mutated()) + if gen >= lastChange + 1000: + print("\n1000 stable iterations; resetting...") + Profile.clear_cache() + newPool = [] + bestThisRound = None + gen = 0 + while len(newPool) < min(poolSize / 2, successes): a = select(candidates, SELECT_P) b = select(candidates, SELECT_P) @@ -199,8 +228,11 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, newPool.append(Profile.random()) pool = newPool gen += 1 + if genLimit is not None and gen >= genLimit: + break except KeyboardInterrupt: + pool.append(bestEver) pool.sort() with open(fileOut, "w") as f: for profile in pool: @@ -214,14 +246,23 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, help='initial altitude (km; defaults to 0)') parser.add_argument('alt1', metavar='alt1', type=int, help='final altitude (km; defaults to the nearest 5 km mark above the atmosphere)') - parser.add_argument('-a', metavar='accel', type=float, default=2, + parser.add_argument('-a', metavar='accel', type=float, default=2.2, help='ship acceleration as a multiple of planet surface gravity (default: %(default)s)') parser.add_argument('-d', metavar='drag', type=float, default=0.2, help='drag coefficient (default: %(default)s)') parser.add_argument('-p', metavar='poolSize', type=int, default=20, help='pool size (default: %(default)s)') + parser.add_argument('--profile', metavar='filename', type=str, + help='profile one generation of execution and save results') #parser.print_help() args = parser.parse_args(sys.argv[1:]) - learnAscent(args.planet, args.alt0, args.alt1, args.a, args.d, args.p) + random.seed(0) + if args.profile: + import cProfile + SILENT = True + cProfile.run('learnAscent(args.planet, args.alt0, args.alt1, args.a, args.d, args.p, genLimit = 5)', args.profile) + else: + SILENT = False + learnAscent(args.planet, args.alt0, args.alt1, args.a, args.d, args.p) From 4ebc842aabaad66b82840efc36158414d7154ba6 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Thu, 28 Feb 2013 20:48:00 -0600 Subject: [PATCH 11/25] Change the default pool size from 20 to 2. --- learn-ascent.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/learn-ascent.py b/learn-ascent.py index 7727ea7..ff0c4f0 100644 --- a/learn-ascent.py +++ b/learn-ascent.py @@ -250,7 +250,7 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, help='ship acceleration as a multiple of planet surface gravity (default: %(default)s)') parser.add_argument('-d', metavar='drag', type=float, default=0.2, help='drag coefficient (default: %(default)s)') - parser.add_argument('-p', metavar='poolSize', type=int, default=20, + parser.add_argument('-p', metavar='poolSize', type=int, default=2, help='pool size (default: %(default)s)') parser.add_argument('--profile', metavar='filename', type=str, help='profile one generation of execution and save results') From 513a21cba22b8293dad312d98201c1742c16f4fc Mon Sep 17 00:00:00 2001 From: gmcnew Date: Sat, 2 Mar 2013 17:16:27 -0600 Subject: [PATCH 12/25] Use Runge-Kutta 4 for better precision. --- ascent.py | 44 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 42 insertions(+), 2 deletions(-) diff --git a/ascent.py b/ascent.py index 921e158..b893f7c 100644 --- a/ascent.py +++ b/ascent.py @@ -239,6 +239,15 @@ def __str__(self): ) v_nothrust = [ v[i] + loss[i] * timestep for i in range(2) ] + def total_accel(pos, vel, thr): + altitude = L2(pos) - planet.radius + ag = planet.gravity(altitude) + ad = planet.drag(altitude, L2(vel), dragCoefficient) + return ( + thr * cos(phi) - (ad * cos(theta) + ag * cos(psi)), + thr * sin(phi) - (ad * sin(theta) + ag * sin(psi)) + ) + def thrustResult2(thrust, tx, ty): return (v_nothrust[0] + thrust * tx, v_nothrust[1] + thrust * ty) @@ -369,10 +378,41 @@ def getAcceleration(dv, alt): thrustLimit = (thrust - a_thrust) if calculateExtraThrust else True thrust = a_thrust + def vmult(val, vector): + return [val * x for x in vector] + + def vadd(v1, v2): + return [v1[i] + v2[i] for i in range(len(v1))] + + def update_rk4(p, v): + dv1 = vmult(timestep, total_accel(p, v, thrust)) + dx1 = vmult(timestep, v) + + dv2 = vmult(timestep, total_accel(vadd(p, vmult(0.5, dx1)), vadd(v, vmult(0.5, dv1)), thrust)) + dx2 = vmult(timestep, vadd(v, vmult(0.5, dv1))) + + dv3 = vmult(timestep, total_accel(vadd(p, vmult(0.5, dx2)), vadd(v, vmult(0.5, dv2)), thrust)) + dx3 = vmult(timestep, vadd(v, vmult(0.5, dv2))) + + dv4 = vmult(timestep, total_accel(vadd(p, dx3), vadd(v, dv3), thrust)) + dx4 = vmult(timestep, vadd(v, dv3)) + + dx = vmult(1.0 / 6.0, vadd(vadd(vadd(dx1, vmult(2, dx2)), vmult(2, dx3)), dx4)) + dv = vmult(1.0 / 6.0, vadd(vadd(vadd(dv1, vmult(2, dv2)), vmult(2, dv3)), dv4)) + + p = vadd(p, dx) + v = vadd(v, dv) + return (p, v) + + def update_euler(p, v): + for i in (0,1): p[i] += v[i] * timestep + v = thrustResult(thrust) + return (p, v) + # Update everything! - for i in (0,1): p[i] += v[i] * timestep + (p, v) = update_rk4(p, v) h = L2(p) - v = thrustResult(thrust) + dV += thrust * timestep dragLoss += a_drag * timestep t += timestep From 438c0697fd388a593d6fe9db16c94d44c84a27bb Mon Sep 17 00:00:00 2001 From: gmcnew Date: Sat, 2 Mar 2013 17:32:24 -0600 Subject: [PATCH 13/25] Fix drag() for negative altitudes on the Mun. --- planet.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/planet.py b/planet.py index 457e5c3..8a2273d 100644 --- a/planet.py +++ b/planet.py @@ -300,7 +300,7 @@ def topOfAtmosphere(self): """ # self.altitude with pressure = self.datumPressure * 1e-6, inlined: # log(1e6) ~ 13.81551... - return self.scale * 13.81551 + return self.scale * 13.81551 if self.scale else -self.radius # the sun has an infinite SOI; arbitrarily set it to 500x the orbit of Jool From eb58aed7aec8aaf3fec3240a9da46f9c5542d782 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Sat, 2 Mar 2013 17:47:54 -0600 Subject: [PATCH 14/25] Add an endAngle parameter to an ascent. --- ascent.py | 8 ++++++-- learn-ascent.py | 51 +++++++++++++++++++++++++++---------------------- 2 files changed, 34 insertions(+), 25 deletions(-) diff --git a/ascent.py b/ascent.py index b893f7c..c4eb736 100644 --- a/ascent.py +++ b/ascent.py @@ -68,7 +68,8 @@ def __init__(self, calculateExtraThrust = False, dragCoefficient = None, specificImpulse = None, - shipThrust = None): + shipThrust = None, + endAngleDeg = 0): """ Compute a climb slope for exiting the atmosphere and achieving orbit. @@ -153,6 +154,9 @@ def __str__(self): if dragCoefficient is None: dragCoefficient = planet.defaultDragCoefficient + # These angles are both measured from horizontal. + endAngleRad = endAngleDeg * math.pi / 180 + # p and v are the current position and velocity; we update these. # We use a coordinate system centered at the core of the planet, # aligned so that the initial altitude forms the y axis. @@ -220,7 +224,7 @@ def __str__(self): ratio **= gravityTurnCurve # The more we did, the closer we want to be to pi/2 from # the straight-up orientation. - phiSurf = - (ratio * math.pi / 2) + phiSurf = - (ratio * (math.pi / 2 - endAngleRad)) # phi is the angle of thrust in the original coordinate frame. phi = phiSurf + psi diff --git a/learn-ascent.py b/learn-ascent.py index ff0c4f0..0a0c48b 100644 --- a/learn-ascent.py +++ b/learn-ascent.py @@ -19,20 +19,22 @@ class Profile: cache = {} MAX_CACHE_SIZE = 10000 - def __init__(self, gt0, gt1, curve): - # Limit values to 3 decimal places of precision. - (gt0, gt1, curve) = [(int(x * 1000) / 1000.0) for x in (gt0, gt1, curve)] + def __init__(self, gt0, gt1, curve, endAngle): + # Limit precision of values. + (gt0, gt1, endAngle) = [(int(x * 100) / 100.0) for x in (gt0, gt1, endAngle)] + curve = int(curve * 1000) / 1000.0 - self.gt0 = sorted([0, gt0, self.alt1])[1] - self.gt1 = sorted([0, gt1, self.alt1])[1] - self.curve = sorted([0, curve, 1])[1] + self.gt0 = sorted([ 0, gt0, self.alt1])[1] + self.gt1 = sorted([ 0, gt1, self.alt1])[1] + self.curve = sorted([ 0, curve, 1])[1] + self.endAngle = sorted([-10, endAngle, 90])[1] - self.score = self.cache.get((self.gt0, self.gt1, self.curve), None) + self.score = self.cache.get((self.gt0, self.gt1, self.curve, self.endAngle), None) if self.score is None: if len(self.cache) == self.MAX_CACHE_SIZE: del self.cache[random.choice(self.cache.keys())] self._calc_score() - self.cache[(self.gt0, self.gt1, self.curve)] = self.score + self.cache[(self.gt0, self.gt1, self.curve, self.endAngle)] = self.score self.generation = 1 @@ -51,10 +53,10 @@ def clear_cache(cls): @classmethod def from_string(cls, str): tokens = [float(f) for f in str.strip().split(" ")] - if len(tokens) == 4: + if len(tokens) == 5: tokens.append(-1) - (gt0, gt1, curve, score) = tokens[:4] - return Profile(gt0, gt1, curve) + (gt0, gt1, curve, endAngle, score) = tokens[:5] + return Profile(gt0, gt1, curve, endAngle) #9.72 23.87 0.3381 1119.067026 @classmethod @@ -62,13 +64,14 @@ def random(cls): gt0 = random.random() * (cls.alt1 / 2) gt1 = gt0 + random.random() * (cls.alt1 - gt0) curve = random.random() * 2 - return Profile(gt0, gt1, curve) + endAngle = random.random() * 90 + return Profile(gt0, gt1, curve, endAngle) def mutated(self): - vals = [self.gt0, self.gt1, self.curve] - i = random.randint(0, 2) + vals = [self.gt0, self.gt1, self.curve, self.endAngle] + i = random.randint(0, len(vals) - 1) vals[i] = self._mutate(vals[i]) - return Profile(vals[0], vals[1], vals[2]) + return Profile(vals[0], vals[1], vals[2], vals[3]) def _mutate(self, value): amount = 0.1 @@ -76,7 +79,7 @@ def _mutate(self, value): # Mutate more slowly in later generations. #amount /= math.log(self.generation + 1) - return value + (random.random() - 0.5) * math.sqrt(value) * amount + return value + (random.random() - 0.5) * (math.sqrt(math.fabs(value)) if value else 1) * amount def _calc_score(self): try: @@ -87,7 +90,8 @@ def _calc_score(self): gravityTurnCurve = self.curve, acceleration = self.planet.gravity() * self.accel, initialAltitude = self.alt0 * 1000, - dragCoefficient = self.drag + dragCoefficient = self.drag, + endAngleDeg = self.endAngle ) self.score = self.ascent.deltaV() except ascent.BadFlightPlanException as bfpe: @@ -102,9 +106,10 @@ def __lt__(self, other): return self.score != -1 and self.score < other.score def combine(self, other): - return Profile(self._combine(self.gt0, other.gt0), - self._combine(self.gt1, other.gt1), - self._combine(self.curve, other.curve)) + return Profile(self._combine(self.gt0, other.gt0), + self._combine(self.gt1, other.gt1), + self._combine(self.curve, other.curve), + self._combine(self.endAngle, other.endAngle)) def better_than(self, other): return (not other) or other.score < 0 or (self.score > 0 and self.score < other.score) @@ -113,7 +118,7 @@ def worse_than(self, other): return (not other) or self.score > other.score def __str__(self): - return "%.3f %.3f %.3f %f" % (self.gt0, self.gt1, self.curve, self.score) + return "%.2f %.2f %.3f %.2f %f" % (self.gt0, self.gt1, self.curve, self.endAngle, self.score) def guide(self): angleStep = 15 @@ -212,8 +217,8 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, newPool.append(candidates[0]) newPool.append(candidates[0].mutated()) - if gen >= lastChange + 1000: - print("\n1000 stable iterations; resetting...") + if gen >= lastChange + 10000: + print("\n10000 stable iterations; resetting...") Profile.clear_cache() newPool = [] bestThisRound = None From b9ea11cecf1a4e858ca64d76e0dfef795b600e55 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Sun, 3 Mar 2013 14:47:17 -0600 Subject: [PATCH 15/25] Clean up ArgumentParser initialization. --- learn-ascent.py | 28 ++++++++++++---------------- 1 file changed, 12 insertions(+), 16 deletions(-) diff --git a/learn-ascent.py b/learn-ascent.py index 0a0c48b..c92c082 100644 --- a/learn-ascent.py +++ b/learn-ascent.py @@ -245,26 +245,22 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, if __name__ == "__main__": parser = argparse.ArgumentParser(description = "Learn an ascent.") - parser.add_argument('planet', metavar='planet', - help='planet or moon') - parser.add_argument('alt0', metavar='alt0', type=int, - help='initial altitude (km; defaults to 0)') - parser.add_argument('alt1', metavar='alt1', type=int, - help='final altitude (km; defaults to the nearest 5 km mark above the atmosphere)') - parser.add_argument('-a', metavar='accel', type=float, default=2.2, - help='ship acceleration as a multiple of planet surface gravity (default: %(default)s)') - parser.add_argument('-d', metavar='drag', type=float, default=0.2, - help='drag coefficient (default: %(default)s)') - parser.add_argument('-p', metavar='poolSize', type=int, default=2, - help='pool size (default: %(default)s)') - parser.add_argument('--profile', metavar='filename', type=str, - help='profile one generation of execution and save results') - #parser.print_help() + args = [ + ('planet', 'planet', str, None, 'planet or moon'), + ('alt0', 'alt0', int, None, 'initial altitude (km)'), + ('alt1', 'alt1', int, None, 'final altitude (km)'), + ('-a', 'accel', float, 2.2, 'ship acceleration as a multiple of planet surface gravity (default: %(default)s)'), + ('-d', 'drag', float, 0.2, 'drag coefficient (default: %(default)s)'), + ('-p', 'poolSize', int, 2, 'pool size (default: %(default)s)'), + ('--profile', 'filename', str, None, 'profile one generation of execution and save results'), + ] + for (name, metavar, type, default, help) in args: + parser.add_argument(name, metavar=metavar, type=type, help=help, default=default) args = parser.parse_args(sys.argv[1:]) - random.seed(0) if args.profile: + random.seed(0) import cProfile SILENT = True cProfile.run('learnAscent(args.planet, args.alt0, args.alt1, args.a, args.d, args.p, genLimit = 5)', args.profile) From 7df47fc6e3bef8bf5fcc2461d75cf3538f75487b Mon Sep 17 00:00:00 2001 From: gmcnew Date: Sun, 3 Mar 2013 18:31:27 -0600 Subject: [PATCH 16/25] Cleanup. --- learn-ascent.py | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/learn-ascent.py b/learn-ascent.py index c92c082..3e715ea 100644 --- a/learn-ascent.py +++ b/learn-ascent.py @@ -7,6 +7,8 @@ import ascent import planet +STABLE_ITERATIONS = 10000 + class Profile: # Class variables. @@ -21,8 +23,8 @@ class Profile: def __init__(self, gt0, gt1, curve, endAngle): # Limit precision of values. - (gt0, gt1, endAngle) = [(int(x * 100) / 100.0) for x in (gt0, gt1, endAngle)] - curve = int(curve * 1000) / 1000.0 + (gt0, gt1, endAngle) = [round(x, 2) for x in (gt0, gt1, endAngle)] + curve = round(curve, 3) self.gt0 = sorted([ 0, gt0, self.alt1])[1] self.gt1 = sorted([ 0, gt1, self.alt1])[1] @@ -161,7 +163,9 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, if fileIn and os.path.exists(fileIn) and not SILENT: with open(fileIn) as f: for line in f.readlines(): - pool.append(Profile.from_string(line)) + line = line.strip() + if line and not line.startswith("#"): + pool.append(Profile.from_string(line)) while len(pool) < poolSize: profile = Profile.random() @@ -217,13 +221,6 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, newPool.append(candidates[0]) newPool.append(candidates[0].mutated()) - if gen >= lastChange + 10000: - print("\n10000 stable iterations; resetting...") - Profile.clear_cache() - newPool = [] - bestThisRound = None - gen = 0 - while len(newPool) < min(poolSize / 2, successes): a = select(candidates, SELECT_P) b = select(candidates, SELECT_P) @@ -232,6 +229,14 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, while len(newPool) < poolSize: newPool.append(Profile.random()) pool = newPool + + if gen >= lastChange + STABLE_ITERATIONS: + print("\n%d stable iterations; resetting..." % STABLE_ITERATIONS) + Profile.clear_cache() + pool = [] + bestThisRound = None + gen = 0 + gen += 1 if genLimit is not None and gen >= genLimit: break From 0228e029f6d61ab6e11b91a8f3282d4298b90391 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Sun, 3 Mar 2013 18:47:15 -0600 Subject: [PATCH 17/25] Cleanup. --- learn-ascent.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/learn-ascent.py b/learn-ascent.py index 3e715ea..79ca019 100644 --- a/learn-ascent.py +++ b/learn-ascent.py @@ -226,17 +226,17 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, b = select(candidates, SELECT_P) newPool.append(a.combine(b)) - while len(newPool) < poolSize: - newPool.append(Profile.random()) - pool = newPool - if gen >= lastChange + STABLE_ITERATIONS: print("\n%d stable iterations; resetting..." % STABLE_ITERATIONS) Profile.clear_cache() - pool = [] + newPool = [] bestThisRound = None gen = 0 + while len(newPool) < poolSize: + newPool.append(Profile.random()) + pool = newPool + gen += 1 if genLimit is not None and gen >= genLimit: break From 0a366afd4c78ec4fff2e0aee94f538e7827c890d Mon Sep 17 00:00:00 2001 From: gmcnew Date: Mon, 25 Mar 2013 20:28:38 -0500 Subject: [PATCH 18/25] Cleanup. --- learn-ascent.py | 8 ++++---- planet.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/learn-ascent.py b/learn-ascent.py index 79ca019..e7cc84d 100644 --- a/learn-ascent.py +++ b/learn-ascent.py @@ -26,10 +26,10 @@ def __init__(self, gt0, gt1, curve, endAngle): (gt0, gt1, endAngle) = [round(x, 2) for x in (gt0, gt1, endAngle)] curve = round(curve, 3) - self.gt0 = sorted([ 0, gt0, self.alt1])[1] - self.gt1 = sorted([ 0, gt1, self.alt1])[1] - self.curve = sorted([ 0, curve, 1])[1] - self.endAngle = sorted([-10, endAngle, 90])[1] + self.gt0 = min(max(gt0, 0), self.alt1) + self.gt1 = min(max(gt1, 0), self.alt1) + self.curve = min(max(curve, 0), 1) + self.endAngle = min(max(endAngle, -10), 90) self.score = self.cache.get((self.gt0, self.gt1, self.curve, self.endAngle), None) if self.score is None: diff --git a/planet.py b/planet.py index 8a2273d..fbf3c5a 100644 --- a/planet.py +++ b/planet.py @@ -335,5 +335,5 @@ def getPlanet(name): return planets[name.lower()] # register the planets with the module, so users can write planet.kerbin or planet.getPlanet("kerbin") equivalently -for (name, p) in planets.iteritems(): +for (name, p) in planets.items(): globals()[name] = p From 73d35c829da0b02184b8e1fddb47bf0efd7e3026 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Mon, 25 Mar 2013 20:30:56 -0500 Subject: [PATCH 19/25] Add recircularize(), which returns delta-V in m/s. --- planet.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/planet.py b/planet.py index fbf3c5a..5874b36 100644 --- a/planet.py +++ b/planet.py @@ -63,6 +63,14 @@ def gravity(self, altitude = 0): r = self.radius + altitude return self.mu / (r*r) + def recircularize(self, alt1, alt2): + (alt1, alt2) = sorted([alt1, alt2]) + v1 = self.orbitalVelocity(alt1) + v2 = self.orbitalVelocity(alt1, ap = alt2) + v3 = self.orbitalVelocity(alt2, pe = alt1) + v4 = self.orbitalVelocity(alt2) + return (v2 - v1) + (v4 - v3) + def orbitalVelocity(self, altitude, ap = None, pe = None): """ Return the velocity required to maintain an orbit with the given From 2626f19eb157703924fc0715f5d3f8510184c413 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Mon, 25 Mar 2013 22:10:54 -0500 Subject: [PATCH 20/25] Assume endAngle should always be 0. --- learn-ascent.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/learn-ascent.py b/learn-ascent.py index e7cc84d..879ed0b 100644 --- a/learn-ascent.py +++ b/learn-ascent.py @@ -31,6 +31,11 @@ def __init__(self, gt0, gt1, curve, endAngle): self.curve = min(max(curve, 0), 1) self.endAngle = min(max(endAngle, -10), 90) + # Calculating the end angle is actually not very helpful. + self.endAngle = 0 + + self.ascent = None + self.score = self.cache.get((self.gt0, self.gt1, self.curve, self.endAngle), None) if self.score is None: if len(self.cache) == self.MAX_CACHE_SIZE: From a2012cd0b4b8c4c545b29e5dcf20086283bd52bc Mon Sep 17 00:00:00 2001 From: gmcnew Date: Mon, 25 Mar 2013 22:19:00 -0500 Subject: [PATCH 21/25] Calculate steering, drag, and gravity losses. --- ascent.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/ascent.py b/ascent.py index c4eb736..b209806 100644 --- a/ascent.py +++ b/ascent.py @@ -193,6 +193,9 @@ def __str__(self): # decide how much to burn and update our position. targetApoapsis = orbitAltitude + planet.radius climbSlope = [] # ClimbPoint list + loss_steering = 0 + loss_drag = 0 + loss_gravity = 0 # Keep track of approximately how much thrust would be needed to reach # the desired apoapsis. @@ -382,6 +385,12 @@ def getAcceleration(dv, alt): thrustLimit = (thrust - a_thrust) if calculateExtraThrust else True thrust = a_thrust + vmag = math.sqrt(v[0] ** 2 + v[1] ** 2) + vdot = 0 if not vmag else (v[0] * cos(phi) + v[1] * sin(phi)) / vmag + loss_steering += timestep * thrust * (1 - vdot) + loss_drag += timestep * a_drag + loss_gravity += timestep * math.fabs(sin(thetaSurf)) * planet.gravity(alt) + def vmult(val, vector): return [val * x for x in vector] @@ -426,6 +435,9 @@ def update_euler(p, v): # Timed out... raise BadFlightPlanException + self.loss_steering = loss_steering + self.loss_drag = loss_drag + self.loss_gravity = loss_gravity self._climbSlope = climbSlope self.planet = planet self.orbit = orbitAltitude From b35f5eff0a54f9eb59931e3fed05e07ff8acec69 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Mon, 25 Mar 2013 22:21:35 -0500 Subject: [PATCH 22/25] Tweak output. --- learn-ascent.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/learn-ascent.py b/learn-ascent.py index 879ed0b..e8b459b 100644 --- a/learn-ascent.py +++ b/learn-ascent.py @@ -129,13 +129,19 @@ def __str__(self): def guide(self): angleStep = 15 - guide = "(angle alt) " - for angle in range(0, 90 + 1, angleStep): - dy = self.gt1 - self.gt0 - alt = self.gt0 + ((float(angle) / 90) ** (1.0 / self.curve)) * dy - if angle != 0: - guide += ", " - guide += "%d %4.1f" % (angle, alt) + guide = "" + """ + if self.curve: + guide = "(angle alt) " + for angle in range(0, 90 + 1, angleStep): + dy = self.gt1 - self.gt0 + alt = self.gt0 + ((float(angle) / 90) ** (1.0 / self.curve)) * dy + if angle != 0: + guide += ", " + guide += "%d %4.1f" % (angle, alt) + """ + if self.ascent: + guide += " loss g/atm/steer=%.2f/%.2f/%.2f" % (self.ascent.loss_gravity, self.ascent.loss_drag, self.ascent.loss_steering) guide += "\n" return guide From 00a7f31cf6f399ca587767196b52a49539d0a0b1 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Tue, 26 Mar 2013 19:50:54 -0500 Subject: [PATCH 23/25] Add impactVelocity() given a free-fall altitude. --- planet.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/planet.py b/planet.py index 5874b36..3700f52 100644 --- a/planet.py +++ b/planet.py @@ -71,6 +71,9 @@ def recircularize(self, alt1, alt2): v4 = self.orbitalVelocity(alt2) return (v2 - v1) + (v4 - v3) + def impactVelocity(self, altitude): + return self.orbitalVelocity(0, ap = altitude, pe = -self.radius) + def orbitalVelocity(self, altitude, ap = None, pe = None): """ Return the velocity required to maintain an orbit with the given From 406c97450e7f937a774d3def1d5cb43cbcbeeeb8 Mon Sep 17 00:00:00 2001 From: gmcnew Date: Fri, 20 Dec 2013 03:07:29 +0000 Subject: [PATCH 24/25] Cleanup. --- learn-ascent.py | 60 ++++++++++++++++++++++++++++++------------------- 1 file changed, 37 insertions(+), 23 deletions(-) diff --git a/learn-ascent.py b/learn-ascent.py index e8b459b..9500a7d 100644 --- a/learn-ascent.py +++ b/learn-ascent.py @@ -7,7 +7,10 @@ import ascent import planet -STABLE_ITERATIONS = 10000 +STABLE_ITERATIONS = 1000 + +# Using an end angle other than 0 is actually not very helpful. +VARY_END_ANGLE = False class Profile: @@ -29,10 +32,7 @@ def __init__(self, gt0, gt1, curve, endAngle): self.gt0 = min(max(gt0, 0), self.alt1) self.gt1 = min(max(gt1, 0), self.alt1) self.curve = min(max(curve, 0), 1) - self.endAngle = min(max(endAngle, -10), 90) - - # Calculating the end angle is actually not very helpful. - self.endAngle = 0 + self.endAngle = min(max(endAngle, -10), 90) if VARY_END_ANGLE else 0 self.ascent = None @@ -140,11 +140,26 @@ def guide(self): guide += ", " guide += "%d %4.1f" % (angle, alt) """ - if self.ascent: - guide += " loss g/atm/steer=%.2f/%.2f/%.2f" % (self.ascent.loss_gravity, self.ascent.loss_drag, self.ascent.loss_steering) - guide += "\n" return guide + @classmethod + def desc_header(cls): + header = "start end shape" + if VARY_END_ANGLE: + header += " endAng" + header += " deltaV" + header += " loss_g atm steer" + return header + + def desc(self): + desc = "" + desc += "%5.2f %6.2f %5.1f " % (self.gt0, self.gt1, self.curve * 100) + if VARY_END_ANGLE: + desc += "%6.2f " % self.endAngle + desc += "%8.2f " % self.score + desc += "%7.2f %7.2f %6.2f" % (self.ascent.loss_gravity, self.ascent.loss_drag, self.ascent.loss_steering) + return desc + def select(pool, s): for profile in pool: if random.random() > s: @@ -161,7 +176,10 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, endAlt = math.ceil(p.topOfAtmosphere() / 5000) * 5 if not SILENT: - print("ascending to %d km with %.2f g's of acceleration" % (endAlt, accel)) + print("ascending on %s from %d to %d km" % (p, startAlt, endAlt)) + print("max acceleration: %.2f x surface gravity = %.2f m/s^2" % (accel, accel * p.gravity())) + if drag != 0.2: + print("drag coefficient: %.2f" % drag) Profile.init(p, startAlt, endAlt, accel, drag) fileOut = "%s_%d_%d_%.2f_%.2f.txt" % (p.name, startAlt, endAlt, accel, drag) @@ -188,29 +206,24 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, needNewline = False bestThisRound = None + if not SILENT: + print("%6s %s" % ("iter", Profile.desc_header())) try: while True: best = None worst = None total = 0 successes = 0 - # 7.6 45.0 0.653 = 4435.43 - # 8.7 45.9 0.610 = 4436.00 candidates = [] for profile in pool: profile.generation = gen - if profile.better_than(best) and not SILENT: + if profile.better_than(best): if profile.better_than(bestEver): bestEver = profile best = profile if profile.better_than(bestThisRound): lastChange = gen bestThisRound = profile - if (needNewline): - sys.stdout.write("\n") - needNewline = False - print(" " * 8 + str(profile)) - print(profile.guide()) if profile.worse_than(worst): worst = profile @@ -220,10 +233,10 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, successes += 1 candidates.append(profile) candidates.sort() - if successes: - needNewline = True - if not SILENT: - sys.stdout.write("\r%6d: best %f, average %f, best ever %f" % (gen, best.score, total / successes, bestEver.score)) + if successes and not SILENT: + lineOut = "\r%6d %s" % (gen, bestEver.desc()) + sys.stdout.write(lineOut) + sys.stdout.flush() newPool = [] SELECT_P = 0.5 @@ -238,11 +251,11 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, newPool.append(a.combine(b)) if gen >= lastChange + STABLE_ITERATIONS: - print("\n%d stable iterations; resetting..." % STABLE_ITERATIONS) + #print("\n%6d stable iterations; resetting..." % STABLE_ITERATIONS) Profile.clear_cache() newPool = [] + lastChange = gen bestThisRound = None - gen = 0 while len(newPool) < poolSize: newPool.append(Profile.random()) @@ -253,6 +266,7 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, break except KeyboardInterrupt: + print("") pool.append(bestEver) pool.sort() with open(fileOut, "w") as f: From 6122806d847f52263842251d36cff7733f23a54e Mon Sep 17 00:00:00 2001 From: gmcnew Date: Fri, 20 Dec 2013 03:55:12 +0000 Subject: [PATCH 25/25] Print an angle/altitude guide on exit. --- learn-ascent.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/learn-ascent.py b/learn-ascent.py index 9500a7d..10e0b38 100644 --- a/learn-ascent.py +++ b/learn-ascent.py @@ -130,16 +130,12 @@ def __str__(self): def guide(self): angleStep = 15 guide = "" - """ if self.curve: - guide = "(angle alt) " + guide = "angle altitude" for angle in range(0, 90 + 1, angleStep): dy = self.gt1 - self.gt0 alt = self.gt0 + ((float(angle) / 90) ** (1.0 / self.curve)) * dy - if angle != 0: - guide += ", " - guide += "%d %4.1f" % (angle, alt) - """ + guide += "\n%5d %8.1f" % (angle, alt) return guide @classmethod @@ -272,6 +268,8 @@ def learnAscent(planetName, startAlt = 0, endAlt = None, accel = 2, drag = 0.2, with open(fileOut, "w") as f: for profile in pool: f.write("%s\n" % profile) + print("") + print(bestEver.guide()) if __name__ == "__main__": parser = argparse.ArgumentParser(description = "Learn an ascent.")