This notebook contains some calculations done for this stackexchange answer for the question how effective is speeding?.
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.misc import imread
%matplotlib inline
I wanted to fit the model I found in the answer: $$ \text{ mpg } = \frac{ \alpha }{ v^2 + \frac{\beta }{v} } $$
To do that I need some data. I found this picture on wikimedia, which we will work off of.
# download the image to /tmp and open with scipy
!wget http://upload.wikimedia.org/wikipedia/commons/a/a9/Fuel_economy_vs_speed_1997.png -O /tmp/Fuel_economy_vs_speed_1997.png
pic = imread('/tmp/Fuel_economy_vs_speed_1997.png')
--2014-08-09 12:54:19-- http://upload.wikimedia.org/wikipedia/commons/a/a9/Fuel_economy_vs_speed_1997.png Resolving upload.wikimedia.org (upload.wikimedia.org)... 208.80.154.240, 2620:0:861:ed1a::2:b Connecting to upload.wikimedia.org (upload.wikimedia.org)|208.80.154.240|:80... connected. HTTP request sent, awaiting response... 200 OK Length: 32006 (31K) [image/png] Saving to: `/tmp/Fuel_economy_vs_speed_1997.png' 100%[======================================>] 32,006 89.4K/s in 0.3s 2014-08-09 12:54:21 (89.4 KB/s) - `/tmp/Fuel_economy_vs_speed_1997.png' saved [32006/32006]
# the image we got
plt.imshow(pic,origin='upper');
Looking at the image, I did my "fit" analytically, since I can calculate, given the model that the maximum mpg occurs when the derivative vanishes
$$ d\text{mpg} = - \alpha ( v^2 + \beta/v)^{-2} ( 2 v - \beta/v^2 ) = 0 $$or at $$ v^3 = \beta/2 \implies \beta = 2 v_\text{peak}^3$$
And at this peak velocity, the peak mpg is $$ \text{mpg}_{\text{peak}} = \frac{ \alpha }{ 3 v_{\text{peak}}^2 } $$ which implies $$ \alpha = 3 v_{\text{peak}}^2 \text{mpg}_{\text{peak}} $$
Looking at the data in the image, I will take as my peak speed of 40 mpg and a peak mpg of 35
# defining the model
beta = 2*(40)**3
alpha = 3*(40)**2*35
model = lambda v: alpha/( v**2 + 2*(40)**3/v )
# just checking that the model looks reasonable
vs = np.linspace(1,90,1000)
mpgs = model(vs)
plt.plot( vs, mpgs);
In order to overlap my picture ontop of the image, I need to get the translation between speed on the graph and mpg on the graph and pixels on the image. So for this I just looked at the graph and for instance saw that the 90 mph scale on the image goes from pixel 100 to 420 or so, so that gives me the conversion factor, and the 0 of the mph scale occurs at pixel 100 or so, which gives me an offset
# looking at hte data, get the conversion factors from mph to pixels and mpg to pixels
vscale = ( 420-100) / 90.
mpgscale = ( 400 - 150 ) / 40.
# transform the model into pixel space
vpic = 100 + vscale*vs
mpgpic = 380 - (mpgscale*mpgs)
#Plot poth on top of one another
plt.figure(figsize=(22,16));
plt.imshow(pic, origin='upper');
plt.plot(vpic, mpgpic, linewidth=30, alpha=0.5);
plt.xlim((0,pic.shape[1]));
plt.ylim((pic.shape[0],0));
plt.axis('off');
#plt.savefig('modified.png')
(0.0, 680.0, 480.0, 0.0)
Here I try to compute the costs of driving both per hour and per distance
gasprice = 3.753 # in dollars from wolfram alpha
costperhour = vs / mpgs * gasprice
plt.plot( vs, costperhour, lw=3 )
plt.grid()
plt.xlabel('speed (mph)')
plt.ylabel('cost ($/hour)');
#plt.savefig('costperhour.png');
plt.title("Operating cost per hour for average US vehicle");
# I get the change by taking a gradient
dcost5mph = np.gradient(costperhour, vs[1]-vs[0])*5
plt.plot(vs, dcost5mph, lw=3)
plt.grid()
plt.xlabel('speed (mph)')
plt.ylabel('dcost ($/hour/5 mph)');
#plt.savefig('dcostperhour.png');
plt.title("Change in operating cost for 5 mph difference in speed per hour");
costpermile = gasprice*10./mpgs
plt.plot(vs[100:], costpermile[100:], lw=3)
plt.grid()
plt.xlabel('speed (mph)')
plt.ylabel('$/10 miles');
plt.title("Cost per 10 miles");
#plt.savefig("costpermile.png");
dcostpermile = np.gradient(costpermile, vs[1]-vs[0])*5;
plt.plot(vs[100:], dcostpermile[100:], lw=3)
plt.grid()
plt.xlabel('speed (mph)')
plt.ylabel(r'delta $/10 miles/5 mph difference');
plt.title("Difference in cost per 10 miles for 5 mph change in speed");
#plt.savefig("dcostpermile.png");
block_size = 1./10 # in miles
speed = 30 # mph
#speed = 35
time_between_lights = block_size / speed * 60 * 60 # in seconds
off_time = 50 # seconds
on_time = 50 # seconds
p_on = on_time *1.0/(on_time + off_time)
n = 100
p_on
0.5
from random import random
def next_block(speed):
time_between_lights = block_size / speed * 60 * 60
if random() < p_on:
return time_between_lights
else:
return time_between_lights + random() * off_time
def difftimes(nblocks, speed1, speed2):
times = np.array([ sum( next_block(speed1) for i in xrange(nblocks) ) for i in xrange(10000) ])
times2 = np.array([ sum( next_block(speed2) for i in xrange(nblocks) ) for i in xrange(10000) ])
return times - times2
bins = np.linspace(-300,300,100)
plt.hist( difftimes(5, 45,50), bins, histtype='step', normed=True, alpha=0.6 ,lw=3, label='N=5');
plt.hist( difftimes(15, 45,50), bins, histtype='step', normed=True, alpha=0.6 ,lw=3, label='N=15');
plt.hist( difftimes(50, 45,50), bins, histtype='step', normed=True, alpha=0.6 ,lw=3, label='N=50');
plt.legend()
plt.xlim((-300,300));
# average speed
tau = 30 # seconds
vs = np.linspace(10,90,1000)
p = 0.65
d = 1./10
average_speed = 60* 60 * d / ( d/vs * 60 * 60 + tau/2 * (1-p) )
plt.plot(vs, average_speed, lw=3);
plt.grid()
plt.xlabel('target speed (mph)');
plt.ylabel('average speed (mph)');
#plt.savefig('averagecityspeed.png');
# Gains in time
plt.plot( vs, 5./vs , lw=3, label="5 mph gain");
plt.xlabel('target speed (mph)');
plt.ylabel('fractional time');
sigma = np.sqrt( tau**2 / 12. * ( 1 - p ) * ( 3 * p - 1) )
mu = tau/2. * ( 1 - p )
avgtime = d/vs * 60 * 60 + mu
plt.plot( vs, sigma/avgtime, ls='--', label='N=1' );
plt.plot( vs, sigma/avgtime/np.sqrt(10), ls='--', label='N=10' );
plt.plot( vs, sigma/avgtime/10.0, ls='--', label='N=100' );
plt.yscale('log');
plt.legend(loc=1);
ax = plt.gca()
plt.grid();
plt.tight_layout()
#plt.savefig('gainsversuslights.png');
plt.plot(vs, avgtime);