From dfacde4c4a6886582d0f64035d3c2a0700766a49 Mon Sep 17 00:00:00 2001 From: Joe Anderson Date: Mon, 31 Oct 2011 12:55:08 -0500 Subject: 5-var optimize, Started writing solver() --- throw.py | 93 ++++++++++++++++++++++++++++++++++++++++------------------------ 1 file changed, 59 insertions(+), 34 deletions(-) diff --git a/throw.py b/throw.py index 2f8d36f..e36b151 100644 --- a/throw.py +++ b/throw.py @@ -65,13 +65,14 @@ cos = taylor_cos sin = taylor_sin ### radii: -global r0 -r0 = 1. #meters -- the length of cable in the first full-circle rotation +global r0, r_roll, length +r0 = 1.151467 #meters -- the length of cable in the first full-circle rotation r_anchor = .2 #meters -- the radius of the anchor point -r_winch = .2 #meters -- the max radius of the roll of cable -length = 1000 #meters -- the length of the cable -phi_max = 2 * (length - r0) / r_winch #rad -r = lambda phi: r0 + r_winch*(1 - phi/(2*phi_max))*phi # cable let out +r_roll = .2045575 #meters -- the max radius of the roll of cable +length = 578.63682 #meters -- the length of the cable +phi_max = 2 * (length - r0) / r_roll #rad +r = lambda phi: r0 + r_roll*(1 - phi/(2*phi_max))*phi # cable length let out +r_dot_co = lambda phi: r_roll*(1 - phi/phi_max) #*phi1 r_sat = lambda phi: numpy.sqrt(r(phi)**2 + r_anchor) # radius of satellite ### gravitational acceleration: @@ -100,7 +101,6 @@ M = lambda phi: Mass - m - l*r(phi) r_cm = lambda phi: (m*r(phi) + .5*l*r(phi)**2) / (m + l*r(phi) + M(phi)) v_cm_co = lambda phi: 1 - r_cm(phi)/r(phi) -tan = lambda phi: r_anchor / r(phi) # tan(theta) def throw(phi, time): """The diff eqs for the windup.""" phi, phi1 = phi @@ -108,26 +108,27 @@ def throw(phi, time): # phi' = phi' phi1, # phi'' = - ( tan(phi) - r_winch/r(phi) ) * phi1**2 - + g(phi) * ( sin(phi) - cos(phi)*tan(phi) ) / r(phi) - - quad_co * r(phi) *(v_cm_co(phi) * phi1 )**2 / m + ( (r_sat(phi)*r_anchor/r(phi) - r_dot_co(phi)) * phi1**2 + + g(phi) * ( sin(phi) - cos(phi)*r_anchor/r(phi) ) + - quad_co * ( v_cm_co(phi) * r(phi) * phi1 )**2 / m + ) / r_sat(phi) ] return prime_vars #good final times: 20 185 1874 -def solve_fun(func, init_vars, space=(0, 1874, 10**5)): +def solve_fun(func, init_vars, timespace): time0 = current_time() - time = numpy.linspace(*space) + time = numpy.linspace(*timespace) soln = odeint(func, init_vars, time) print "took %s seconds to solve" % (current_time() - time0) return time, soln - -def interpret(soln, max_tension=max_tension): + +TimespaceTooSmall = Warning(" *** Requires larger timespace!!! *** ") +def interpret(soln, max_tension): phi = soln[:, 0] phi1 = soln[:, 1] v_tan = v_cm_co(phi) * phi1 * r(phi) / 1000 #kilometers/second - cos_t = lambda phi: cos(numpy.arctan(tan(phi))) # cos(theta) - tension = (m/cos_t(phi))*(phi1**2 * r(phi) - g(tau/4))/1000 #kiloNewtons + tension = (m*r_sat(phi)/r(phi))*(phi1**2 * r(phi) - g(tau/4))/1000 #kN cable_moment = l*((r(phi) - r_cm(phi))**3 - r_cm(phi)**3)/3 #kg m^2 moment = m*(r(phi) - r_cm(phi))**2 + cable_moment - M(phi)*r_cm(phi)**2 momentum = moment*phi1**2 #kg m^2/s^2 @@ -137,14 +138,32 @@ def interpret(soln, max_tension=max_tension): max_arg = numpy.argmax(v_tan[:maxindex]) except ValueError: max_arg = numpy.argmax(v_tan) - return phi, phi1, v_tan, cos_t, tension, moment, momentum, max_arg + if phi[max_arg] > abs(phi[-1] - tau): + raise TimespaceTooSmall + + return phi, phi1, v_tan, tension, moment, momentum, max_arg + +def solver(func, init_vars, timespace=(0,600,10**5), max_tension=max_tension): + for n in itertools.count(1): + time, soln = solve_fun(func, init_vars, timespace) + + try: phi, phi1, v_tan, tension, moment, momentum, max_arg = \ + interpret(soln, max_tension) + except TimespaceTooSmall: + timespace = map(lambda x: x*.2*1.5**n, timespace) + continue + break + + phi, phi1, v_tan, tension, moment, momentum, time = \ + map(lambda x: numpy.ndarray.__getslice__(x, 0, max_arg+1), + [phi, phi1, v_tan, tension, moment, momentum, time]) def plot(time, soln, graphs=['all'], i=7): - phi, phi1, v_tan, cos_t, tension, moment, momentum, max_arg = \ + phi, phi1, v_tan, tension, moment, momentum, max_arg = \ interpret(soln) - phi_max, phi1_max, v_tan_max, tension_max, moment_max, momentum_max = \ - map(lambda x: numpy.ndarray.__getitem__(x, max_arg), - [phi, phi1, v_tan, tension, moment, momentum]) + (phi_max, phi1_max, v_tan_max, tension_max, moment_max, momentum_max, + time_max) = map(lambda x: numpy.ndarray.__getitem__(x, max_arg), + [phi, phi1, v_tan, tension, moment, momentum, time]) def test(*args): b = 0 @@ -223,28 +242,33 @@ def plot(time, soln, graphs=['all'], i=7): def opt_fun(guesses, f, func, init_vars): """Solve, plot, and compare the diff eq.""" - global r0, m, Mass - r0, m, Mass = guesses + global r0, r_roll, length, m, Mass + r0, r_roll, length, m, Mass = guesses time, soln = solve_fun(func, init_vars) - phi, phi1, v_tan, cos_t, tension, moment, momentum, max_arg = \ + phi, phi1, v_tan, tension, moment, momentum, max_arg = \ interpret(soln) - phi_max, phi1_max, v_tan_max, tension_max, moment_max, momentum_max = \ - map(lambda x: numpy.ndarray.__getitem__(x, max_arg), - [phi, phi1, v_tan, tension, moment, momentum]) + phi, phi1, v_tan, tension, moment, momentum, time = \ + map(lambda x: numpy.ndarray.__getslice__(x, 0, max_arg+1), + [phi, phi1, v_tan, tension, moment, momentum, time]) + # Parameter values. + f.write("%s, " % r_roll) + f.write("%s, " % length) f.write("%s, " % r0) f.write("%s, " % m) f.write("%s, " % Mass) - f.write("%s, " % phi_max) - f.write("%s, " % phi1_max) - f.write("%s, " % tension_max) - f.write("%s\n" % v_tan_max) + + # Values at the max velocity point. + f.write("%s, " % phi[-1]) + f.write("%s, " % phi1[-1]) + f.write("%s, " % tension[-1]) + f.write("%s\n" % v_tan[-1]) # If you want strict limits, you can change the return to: #return 1/v_tan + 1000*(Mass > Mass_limit) + 1000*(phi1.max() > phi1_limit) - return 2/v_tan_max + Mass/600 + phi1.max() / tau + return 2/v_tan[-1] + Mass/600 + phi1[-1] / tau v0 = numpy.sqrt(g(0) * r(0)) w0 = v0/r0 @@ -254,8 +278,9 @@ inits = [0, w0] outfile = '/home/jandew/Documents/Impossible Challenges/microsatellite/out' def optimize(): fil = file(outfile, 'w') - fil.write('r0, sat. mass, Total Mass, phi, omega, Tension, Max Vel.\n') - guess = (r0, m, Mass) + fil.write( + 'r0, r_roll, length of cable, m, Mass, phi, omega, Tension, Max Vel.\n') + guess = (r0, r_roll, length, m, Mass) #try: print opt_fun(guess, fil, throw, inits) try: print fmin(opt_fun, guess, args=(fil, throw, inits)) finally: fil.close() -- cgit v1.2.3-54-g00ecf