Using GPkit for model representation and solver interfacing.
An airplane with that consumes as little fuel as possible when cruising (while remaining capable of taking off from the ground, with wings that won't break, etc.)
Start by importing Numpy and GPkit, and turning on LaTex printing for GPkit variables and equations:
In [1]:
import numpy as np
import gpkit
gpkit.init_ipynb_printing()
Then define the Hermite interpolation function for 1-D Hermite splines:
In [2]:
def hermite_interp(x, y, m, xvec):
'''
Evaluate the piecewise cubic Hermite interpolant with monoticity preserved
x = array containing the x-data
y = array containing the y-data
m = slopes at each (x,y) point [computed to preserve monotonicity]
xnew = new "x" value where the interpolation is desired
x must be sorted low to high... (no repeats)
y can have repeated values
This works with either a scalar or vector of "xvec"
'''
n = len(x)
mm = len(xvec)
# Create "copies" of "x" as rows in a mxn 2-dimensional vector
xx = np.resize(x,(mm,n)).transpose()
xxx = xx > xvec
# Compute column by column differences
z = xxx[:-1,:] - xxx[1:,:]
# Collapse over rows...
k = z.argmax(axis=0)
# Create the Hermite coefficients
h = x[k+1] - x[k]
t = (xvec - x[k]) / h
# Hermite basis functions
h00 = (2 * t**3) - (3 * t**2) + 1
h10 = t**3 - (2 * t**2) + t
h01 = (-2* t**3) + (3 * t**2)
h11 = t**3 - t**2
# Compute the interpolated value of "y"
ynew = h00*y[k] + h10*h*m[k] + h01*y[k+1] + h11*h*m[k+1]
return ynew
In [3]:
# For scipyp n-dimensional interpolation (not working yet)
from scipy.interpolate import interpnd
Declare variables,
In [4]:
free_variables = {
'A': "aspect ratio",
'S': ["m^2", "wing area"],
'C_D': "wing drag coefficient",
'C_L': "wing lift coefficent",
'C_f': "skin friction coefficient",
'Re': "Reynold's number",
'W': ["N", "aircraft weight"],
'W_w': ["N", "wing weight"],
'D': ["N", "cruising drag"],
'V': ["m/s", "cruising speed"],
'V_min': ["m/s", "takeoff speed"],
}
gpkit.monify_up(globals(), free_variables)
and constants:
In [5]:
constants = {
'pi': (np.pi, "half of the circle constant"),
'CDA0': (0.03062702, "m^2", "fuselage drag area"),
'rho': (1.23, "kg m^3", "density of air"),
'mu': (1.78e-5, "kg/m*s", "viscosity of air"),
'S_wetratio': (2.05, "wetted area ratio"),
'k': (1.2, "form factor"),
'e': (0.96, "Oswald efficiency factor"),
'W_0': (4940, "N", "aircraft weight excluding wing"),
'N_ult': (2.5, "ultimate load factor"),
'tau': (0.12, "airfoil thickness to chord ratio"),
'C_Lmax': (2.0, "max CL with flaps down"),
}
gpkit.monify_up(globals(), constants)
In [6]:
equations = []
C_D_fuse = CDA0/S # Fuselage viscous drag,
C_D_wpar = k*C_f*S_wetratio # wing parasitic drag,
C_D_ind = C_L**2/(pi*A*e) # and induced drag
equations += [C_D == C_D_fuse + C_D_wpar + C_D_ind] # together make a model for overall drag
W_w_strc = 8.71e-5*(N_ult*A**1.5*(W_0*W*S)**0.5)/tau # A model for the wing's structural members
W_w_surf = 45.24*S # and the weight of the wing's 'skin'
equations += [W_w == W_w_surf + W_w_strc] # form a model of total wing weight.
Then create the whole GP:
In [7]:
gp = gpkit.GP( # Minimize
D, # cruising drag
[ # subject to
D >= 0.5*rho*S*C_D*V**2, # the definition of drag,
Re <= (rho/mu)*V*(S/A)**0.5, # the definition of Reynold's number,
C_f >= 0.074/Re**0.2, # a skin friction model,
W <= 0.5*rho*S*C_L*V**2, # a cruising lift constraint,
W <= 0.5*rho*S*C_Lmax*V_min**2, # a takeoff lift constraint,
W == W_0 + W_w, # the definition of aircraft weight,
] + equations, substitutions=constants) # and the equations and constants we defined earlier.
In [8]:
gp.presweep = gp.last
Load plotting tools:
In [9]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams["figure.frameon"] = False
matplotlib.rcParams["legend.frameon"] = False
Declare some helper functions:
In [10]:
def dict_append(base, new):
if not base:
base.update(new)
for key, item in base.items():
if type(item) is not list:
base[key] = [base[key]]
base[key].append(new[key])
def _autosweep(xs, cs, dc_dxs, data):
"""
Given two positions and slopes, guesses at the intersection and checks.
If the guess was within eps%, end: else, split
"""
global i, n_dims, eps
A = np.array([[-1] + [dc_dxs[row][col]
for col in range(n_dims)]
for row in range(n_dims+1)])
b = np.array([np.inner(dc_dxs[row], xs[row]) - cs[row]
for row in range(n_dims+1)])
s = np.linalg.solve(A,b)
c_i = s[0]
x = s[1:]
# In the 1-D case intersections are guaranteed to be inside the swept area
# but this does not appear to be the case for higher dimensions.
# (shouldn't it be, though?)
if any(x > xs.max(axis=0)) or any(x < xs.min(axis=0)):
x = xs.sum(axis=0)/(1.0 + n_dims)
if n_dims == 1:
c_g = hermite_interp(xs, cs, dc_dxs, x)
elif n_dims == 2:
# Doesn't work yet; the guesses are bad
cti = interpnd.CloughTocher2DInterpolator(xs, cs)
# and using the actual gradient makes them inaccurate -and- negative
cti.grad = np.array([ [dc_dx] for dc_dx in dc_dxs])
c_g = cti(x)
else:
# Untested
cti = interpnd.LinearNDInterpolator(xs, cs)
c_g = cti(x)
gp.load(gp.presweep, print_boundwarnings=False)
gp.sub(dict(zip(autosweep.keys(), x)))
newdata = gp.solve(printing=False)
dict_append(data, newdata)
c_a = newdata['D']
dc_dx = np.array([newdata["S{%s}"%svar] for svar in autosweep]) * newdata["D"]/x
if n_dims == 1:
plt.plot([xs[0], x], [cs[0], c_i], 'k', alpha=0.1, linewidth=2)
plt.plot([xs[1], x], [cs[1], c_i], 'k', alpha=0.1, linewidth=2)
plt.plot(x, c_a, 'ko', alpha=0.5, markeredgecolor='none')
plt.plot(x, c_g, 'ro', alpha=0.5, markeredgecolor='none')
elif n_dims == 2:
plt.plot(x[0], x[1], 'ko', alpha=0.5, linewidth=0)
#print xs, x
#print cs
#print "guessed %.2g, but it was actually %.2g" % (c_g, c_a)
if abs(1-c_a/c_g) > 10*eps**2:
for j in range(n_dims+1):
data = _autosweep(np.array(list(xs[:j])+[x]+list(xs[j+1:])),
np.array(list(cs[:j])+[c_a]+list(cs[j+1:])),
np.array(list(dc_dxs[:j])+[dc_dx]+list(dc_dxs[j+1:])), data)
i += 1
return data
Determine sweep parameters:
In [11]:
autosweep = {"V_min": [10, 50]}
n_dims = len(autosweep)
if n_dims == 1:
sweep_range = autosweep.values()[0]
xs = [[sweep_range[0]], [sweep_range[1]]]
else:
xs = list(np.array(np.meshgrid(*autosweep.values())).T.reshape((2**n_dims, n_dims)))
cs, dc_dxs = [], []
startdata = {}
for x in xs:
gp.load(gp.presweep, print_boundwarnings=False)
gp.sub(dict(zip(autosweep.keys(), x)))
newdata = gp.solve(printing=False)
cs.append(newdata["D"])
dc_dxs.append(np.array([newdata["S{%s}"%svar] for svar in autosweep]) * newdata["D"]/x)
dict_append(startdata, newdata)
xs, cs, dc_dxs = map(np.array, [xs, cs, dc_dxs])
print xs
and ground-truth:
In [12]:
gp.load(gp.presweep, print_boundwarnings=False)
gp.sub({'V_min': ('sweep', np.linspace(10, 50, 100), "m/s", "takeoff speed")})
sweptgp = gp.solve()
gp.sweep = {}
Then autosweep between the given bounds:
In [13]:
import time
eps, i = 0.1, 0
des_errs = 10**(-np.array(range(1, 5), 'f'))
sols = []
act_errs = []
for des_err in des_errs:
start_time = time.time()
eps = des_err
i = n_dims + 1
data = dict(startdata)
plt.figure()
data = _autosweep(xs[:n_dims+1], cs[:n_dims+1], dc_dxs[:n_dims+1], data)
time_taken = time.time()-start_time
x_s, y_s = sweptgp["V_min"], sweptgp["D"]
plt.plot(x_s, y_s, 'k--')
x_as, y_as, S_as = map(np.array, [data["V_min"], data["D"], data["S{V_min}"]])
sort_I = np.argsort(x_as)
x_as, y_as, S_as = x_as[sort_I], y_as[sort_I], S_as[sort_I]
y_asi = hermite_interp(x_as, y_as, S_as*y_as/x_as, x_s[:-1])
plt.plot(x_s[:-1], y_asi, 'r--')
worst_error = abs(y_asi/y_s[:-1] - 1).max()
plt.title("[%i solutions and %.2g seconds] \n Desired error bound was %.0e; we got %.0e"
% (i, time_taken, eps, worst_error))
plt.grid()
plt.ylabel("Drag [N]")
plt.xlabel("V_min [m/s]")
sols.append(i)
act_errs.append(worst_error)
In [14]:
plt.loglog(sols, np.array(des_errs), 'r--')
plt.loglog(sols, (10*np.array(des_errs)**2), 'r')
plt.loglog(sols, np.array(act_errs), 'k')
plt.legend(["Requested Tolerance", "Worst Seen", "Worst Actual"])
plt.xlabel("Inverse Error")
plt.xlabel("Solutions")
lsr = np.log(sols[-1]/sols[0])
print "Log-slope of desired error secant: %.2g" % (np.log(des_errs[-1]/des_errs[0])/lsr)
print " Log-slope of actual error secant: %.2g" % (np.log(act_errs[-1]/act_errs[0])/lsr)
print " Log-slope of seen error secant: %.2g" % ((np.log(10)+2*np.log(act_errs[-1]/act_errs[0]))/lsr)
In the graph above we can also see that the number of solutions per desired error bound follows a one-quarter power law, as might be expected for evenly-spaced spline points [1].
However, actually accuracy lags behind, following a square-root. To correct for this, the desired tolerance is squared and then multiplied by 10: as can be seen, this does a pretty good job of following the actual error.
This example was written by Edward Burnell (eburn@mit.edu), who welcomes any comments or suggestions.