Examples¶
iPython Notebook Examples¶
More examples, including some with in-depth explanations and interactive visualizations, can be seen on nbviewer.
A Trivial GP¶
The most trivial GP we can think of: minimize \(x\) subject to the constraint \(x \ge 1\).
"Very simple problem: minimize x while keeping x greater than 1."
from gpkit import Variable, Model
# Decision variable
x = Variable('x')
# Constraint
constraints = [x >= 1]
# Objective (to minimize)
objective = x
# Formulate the Model
m = Model(objective, constraints)
# Solve the Model
sol = m.solve(verbosity=0)
# print selected results
print("Optimal cost: %s" % sol['cost'])
print("Optimal x val: %s" % sol(x))
Of course, the optimal value is 1. Output:
Optimal cost: 1.0
Optimal x val: 1.0
Maximizing the Volume of a Box¶
This example comes from Section 2.4 of the GP tutorial, by S. Boyd et. al.
"Maximizes box volume given area and aspect ratio constraints."
from gpkit import Variable, Model
# Parameters
alpha = Variable("alpha", 2, "-", "lower limit, wall aspect ratio")
beta = Variable("beta", 10, "-", "upper limit, wall aspect ratio")
gamma = Variable("gamma", 2, "-", "lower limit, floor aspect ratio")
delta = Variable("delta", 10, "-", "upper limit, floor aspect ratio")
A_wall = Variable("A_{wall}", 200, "m^2", "upper limit, wall area")
A_floor = Variable("A_{floor}", 50, "m^2", "upper limit, floor area")
# Decision variables
h = Variable("h", "m", "height")
w = Variable("w", "m", "width")
d = Variable("d", "m", "depth")
#Constraints
constraints = [A_wall >= 2*h*w + 2*h*d,
A_floor >= w*d,
h/w >= alpha,
h/w <= beta,
d/w >= gamma,
d/w <= delta]
#Objective function
V = h*w*d
objective = 1/V # To maximize V, we minimize its reciprocal
# Formulate the Model
m = Model(objective, constraints)
# Solve the Model and print the results table
sol = m.solve(verbosity=0)
print sol.table()
The output is
Cost
----
0.003674 [1/m**3]
Free Variables
--------------
d : 8.17 [m] depth
h : 8.163 [m] height
w : 4.081 [m] width
Constants
---------
A_{floor} : 50 [m**2] upper limit, floor area
A_{wall} : 200 [m**2] upper limit, wall area
alpha : 2 lower limit, wall aspect ratio
beta : 10 upper limit, wall aspect ratio
delta : 10 upper limit, floor aspect ratio
gamma : 2 lower limit, floor aspect ratio
Sensitivities
-------------
alpha : 0.5 lower limit, wall aspect ratio
A_{wall} : -1.5 upper limit, wall area
Water Tank¶
Say we had a fixed mass of water we wanted to contain within a tank, but also wanted to minimize the cost of the material we had to purchase (i.e. the surface area of the tank):
"Minimizes cylindrical tank surface area for a particular volume."
from gpkit import Variable, VectorVariable, Model
M = Variable("M", 100, "kg", "Mass of Water in the Tank")
rho = Variable("\\rho", 1000, "kg/m^3", "Density of Water in the Tank")
A = Variable("A", "m^2", "Surface Area of the Tank")
V = Variable("V", "m^3", "Volume of the Tank")
d = VectorVariable(3, "d", "m", "Dimension Vector")
constraints = (A >= 2*(d[0]*d[1] + d[0]*d[2] + d[1]*d[2]),
V == d[0]*d[1]*d[2],
M == V*rho)
m = Model(A, constraints)
sol = m.solve(verbosity=0)
print sol.table()
The output is
Cost
----
1.293 [m**2]
Free Variables
--------------
A : 1.293 [m**2] Surface Area of the Tank
V : 0.1 [m**3] Volume of the Tank
\vec{d} : [ 0.464 0.464 0.464 ] [m] Dimension Vector
Constants
---------
M : 100 [kg] Mass of Water in the Tank
\rho : 1000 [kg/m**3] Density of Water in the Tank
Sensitivities
-------------
M : 0.6667 Mass of Water in the Tank
\rho : -0.6667 Density of Water in the Tank
Simple Wing¶
This example comes from Section 3 of Geometric Programming for Aircraft Design Optimization, by W. Hoburg and P. Abbeel.
"Minimizes airplane drag for a simple drag and structure model."
import numpy as np
from gpkit import Variable, Model
# Constants
k = Variable("k", 1.2, "-", "form factor")
e = Variable("e", 0.95, "-", "Oswald efficiency factor")
mu = Variable("\\mu", 1.78e-5, "kg/m/s", "viscosity of air")
pi = Variable("\\pi", np.pi, "-", "half of the circle constant")
rho = Variable("\\rho", 1.23, "kg/m^3", "density of air")
tau = Variable("\\tau", 0.12, "-", "airfoil thickness to chord ratio")
N_ult = Variable("N_{ult}", 3.8, "-", "ultimate load factor")
V_min = Variable("V_{min}", 22, "m/s", "takeoff speed")
C_Lmax = Variable("C_{L,max}", 1.5, "-", "max CL with flaps down")
S_wetratio = Variable("(\\frac{S}{S_{wet}})", 2.05, "-", "wetted area ratio")
W_W_coeff1 = Variable("W_{W_{coeff1}}", 8.71e-5, "1/m",
"Wing Weight Coefficent 1")
W_W_coeff2 = Variable("W_{W_{coeff2}}", 45.24, "Pa",
"Wing Weight Coefficent 2")
CDA0 = Variable("(CDA0)", 0.031, "m^2", "fuselage drag area")
W_0 = Variable("W_0", 4940.0, "N", "aircraft weight excluding wing")
# Free Variables
D = Variable("D", "N", "total drag force")
A = Variable("A", "-", "aspect ratio")
S = Variable("S", "m^2", "total wing area")
V = Variable("V", "m/s", "cruising speed")
W = Variable("W", "N", "total aircraft weight")
Re = Variable("Re", "-", "Reynold's number")
C_D = Variable("C_D", "-", "Drag coefficient of wing")
C_L = Variable("C_L", "-", "Lift coefficent of wing")
C_f = Variable("C_f", "-", "skin friction coefficient")
W_w = Variable("W_w", "N", "wing weight")
constraints = []
# Drag model
C_D_fuse = CDA0/S
C_D_wpar = k*C_f*S_wetratio
C_D_ind = C_L**2/(pi*A*e)
constraints += [C_D >= C_D_fuse + C_D_wpar + C_D_ind]
# Wing weight model
W_w_strc = W_W_coeff1*(N_ult*A**1.5*(W_0*W*S)**0.5)/tau
W_w_surf = W_W_coeff2 * S
constraints += [W_w >= W_w_surf + W_w_strc]
# and the rest of the models
constraints += [D >= 0.5*rho*S*C_D*V**2,
Re <= (rho/mu)*V*(S/A)**0.5,
C_f >= 0.074/Re**0.2,
W <= 0.5*rho*S*C_L*V**2,
W <= 0.5*rho*S*C_Lmax*V_min**2,
W >= W_0 + W_w]
print("SINGLE\n======")
m = Model(D, constraints)
sol = m.solve(verbosity=0)
print(sol.table())
print("SWEEP\n=====")
N = 2
sweeps = {V_min: ("sweep", np.linspace(20, 25, N)),
V: ("sweep", np.linspace(45, 55, N)), }
m.substitutions.update(sweeps)
sweepsol = m.solve(verbosity=0)
print(sweepsol.table())
The output is
SINGLE
======
Cost
----
303.1 [N]
Free Variables
--------------
A : 8.46 aspect ratio
C_D : 0.02059 Drag coefficient of wing
C_L : 0.4988 Lift coefficent of wing
C_f : 0.003599 skin friction coefficient
D : 303.1 [N] total drag force
Re : 3.675e+06 Reynold's number
S : 16.44 [m**2] total wing area
V : 38.15 [m/s] cruising speed
W : 7341 [N] total aircraft weight
W_w : 2401 [N] wing weight
Constants
---------
(CDA0) : 0.031 [m**2] fuselage drag area
(\frac{S}{S_{wet}}) : 2.05 wetted area ratio
C_{L,max} : 1.5 max CL with flaps down
N_{ult} : 3.8 ultimate load factor
V_{min} : 22 [m/s] takeoff speed
W_0 : 4940 [N] aircraft weight excluding wing
W_{W_{coeff1}} : 8.71e-05 [1/m] Wing Weight Coefficent 1
W_{W_{coeff2}} : 45.24 [Pa] Wing Weight Coefficent 2
\mu : 1.78e-05 [kg/m/s] viscosity of air
\pi : 3.142 half of the circle constant
\rho : 1.23 [kg/m**3] density of air
\tau : 0.12 airfoil thickness to chord ratio
e : 0.95 Oswald efficiency factor
k : 1.2 form factor
Sensitivities
-------------
W_0 : 1.011 aircraft weight excluding wing
k : 0.4299 form factor
(\frac{S}{S_{wet}}) : 0.4299 wetted area ratio
W_{W_{coeff1}} : 0.2903 Wing Weight Coefficent 1
N_{ult} : 0.2903 ultimate load factor
W_{W_{coeff2}} : 0.1303 Wing Weight Coefficent 2
(CDA0) : 0.09156 fuselage drag area
\mu : 0.08599 viscosity of air
C_{L,max} : -0.1839 max CL with flaps down
\rho : -0.2269 density of air
\tau : -0.2903 airfoil thickness to chord ratio
V_{min} : -0.3678 takeoff speed
e : -0.4785 Oswald efficiency factor
\pi : -0.4785 half of the circle constant
SWEEP
=====
Cost
----
[ 338 294 396 326 ] [N]
Sweep Variables
---------------
V : [ 45 45 55 55 ] [m/s] cruising speed
V_{min} : [ 20 25 20 25 ] [m/s] takeoff speed
Free Variables
--------------
A : [ 6.2 8.84 4.77 7.16 ] aspect ratio
C_D : [ 0.0146 0.0196 0.0123 0.0157 ] Drag coefficient of wing
C_L : [ 0.296 0.463 0.198 0.31 ] Lift coefficent of wing
C_f : [ 0.00333 0.00361 0.00314 0.00342 ] skin friction coefficient
D : [ 338 294 396 326 ] [N] total drag force
Re : [ 5.38e+06 3.63e+06 7.24e+06 4.75e+06 ] Reynold's number
S : [ 18.6 12.1 17.3 11.2 ] [m**2] total wing area
W : [ 6.85e+03 6.97e+03 6.4e+03 6.44e+03 ] [N] total aircraft weight
W_w : [ 1.91e+03 2.03e+03 1.46e+03 1.5e+03 ] [N] wing weight
Constants
---------
(CDA0) : 0.031 [m**2] fuselage drag area
(\frac{S}{S_{wet}}) : 2.05 wetted area ratio
C_{L,max} : 1.5 max CL with flaps down
N_{ult} : 3.8 ultimate load factor
W_0 : 4940 [N] aircraft weight excluding wing
W_{W_{coeff1}} : 8.71e-05 [1/m] Wing Weight Coefficent 1
W_{W_{coeff2}} : 45.24 [Pa] Wing Weight Coefficent 2
\mu : 1.78e-05 [kg/m/s] viscosity of air
\pi : 3.142 half of the circle constant
\rho : 1.23 [kg/m**3] density of air
\tau : 0.12 airfoil thickness to chord ratio
e : 0.95 Oswald efficiency factor
k : 1.2 form factor
Sensitivities
-------------
W_0 : [ 0.919 0.947 0.845 0.847 ] aircraft weight excluding wing
V : [ 0.589 0.249 0.975 0.746 ] cruising speed
k : [ 0.561 0.454 0.63 0.536 ] form factor
(\frac{S}{S_{wet}}) : [ 0.561 0.454 0.63 0.536 ] wetted area ratio
W_{W_{coeff1}} : [ 0.179 0.247 0.108 0.155 ] Wing Weight Coefficent 1
N_{ult} : [ 0.179 0.247 0.108 0.155 ] ultimate load factor
(CDA0) : [ 0.114 0.131 0.146 0.177 ] fuselage drag area
W_{W_{coeff2}} : [ 0.141 0.0911 0.126 0.0787 ] Wing Weight Coefficent 2
\mu : [ 0.112 0.0907 0.126 0.107 ] viscosity of air
\rho : [ -0.172 -0.129 -0.097 -0.0331 ] density of air
\tau : [ -0.179 -0.247 -0.108 -0.155 ] airfoil thickness to chord ratio
e : [ -0.325 -0.415 -0.225 -0.287 ] Oswald efficiency factor
\pi : [ -0.325 -0.415 -0.225 -0.287 ] half of the circle constant
C_{L,max} : [ -0.411 -0.207 -0.521 -0.353 ] max CL with flaps down
V_{min} : [ -0.822 -0.415 -1.04 -0.705 ] takeoff speed
Simple Beam¶
In this example we consider a beam subjected to a uniformly distributed transverse force along its length. The beam has fixed geometry so we are not optimizing its shape, rather we are simply solving a discretization of the Euler-Bernoulli beam bending equations using GP.
"""
A simple beam example with fixed geometry. Solves the discretized
Euler-Bernoulli beam equations for a constant distributed load
"""
import numpy as np
from gpkit import Variable, VectorVariable, Model, ureg
from gpkit.small_scripts import mag
class Beam(Model):
"""Discretization of the Euler beam equations for a distributed load.
Arguments
---------
N : int
Number of finite elements that compose the beam.
L : float
[m] Length of beam.
EI : float
[N m^2] Elastic modulus times cross-section's area moment of inertia.
q : float or N-vector of floats
[N/m] Loading density: can be specified as constants or as an array.
"""
def __init__(self, N=4, **kwargs):
EI = Variable("EI", 1e4, "N*m^2")
dx = Variable("dx", "m", "Length of an element")
L = Variable("L", 5, "m", "Overall beam length")
q = VectorVariable(N, "q", 100*np.ones(N), "N/m",
"Distributed load at each point")
V = VectorVariable(N, "V", "N", "Internal shear")
V_tip = Variable("V_{tip}", 0, "N", "Tip loading")
M = VectorVariable(N, "M", "N*m", "Internal moment")
M_tip = Variable("M_{tip}", 0, "N*m", "Tip moment")
th = VectorVariable(N, "\\theta", "-", "Slope")
th_base = Variable("\\theta_{base}", 0, "-", "Base angle")
w = VectorVariable(N, "w", "m", "Displacement")
w_base = Variable("w_{base}", 0, "m", "Base deflection")
# below: trapezoidal integration to form a piecewise-linear
# approximation of loading, shear, and so on
# shear and moment increase from tip to base (left > right)
shear_eq = (V >= V.right + 0.5*dx*(q + q.right))
shear_eq[-1] = (V[-1] >= V_tip) # tip boundary condition
moment_eq = (M >= M.right + 0.5*dx*(V + V.right))
moment_eq[-1] = (M[-1] >= M_tip)
# slope and displacement increase from base to tip (right > left)
theta_eq = (th >= th.left + 0.5*dx*(M + M.left)/EI)
theta_eq[0] = (th[0] >= th_base) # base boundary condition
displ_eq = (w >= w.left + 0.5*dx*(th + th.left))
displ_eq[0] = (w[0] >= w_base)
# minimize tip displacement (the last w)
Model.__init__(self, w[-1],
[shear_eq, moment_eq, theta_eq, displ_eq,
L == (N-1)*dx], **kwargs)
b = Beam(N=6, substitutions={"L": 6, "EI": 1.1e4, "q": 110*np.ones(6)})
b.zero_lower_unbounded_variables()
sol = b.solve(verbosity=0)
print sol.table()
w_gp = sol("w") # deflection along beam
L, EI, q = sol("L"), sol("EI"), sol("q")
x = np.linspace(0, mag(L), len(q))*ureg.m # position along beam
q = q[0] # assume uniform loading for the check below
w_exact = q/(24.*EI) * x**2 * (x**2 - 4*L*x + 6*L**2) # analytic soln
assert max(abs(w_gp - w_exact)) <= 1.1*ureg.cm
PLOT = False
if PLOT:
import matplotlib.pyplot as plt
x_exact = np.linspace(0, L, 1000)
w_exact = q/(24.*EI) * x_exact**2 * (x_exact**2 - 4*L*x_exact + 6*L**2)
plt.plot(x, w_gp, color='red', linestyle='solid', marker='^',
markersize=8)
plt.plot(x_exact, w_exact, color='blue', linestyle='dashed')
plt.xlabel('x [m]')
plt.ylabel('Deflection [m]')
plt.axis('equal')
plt.legend(['GP solution', 'Analytical solution'])
plt.show()
The output is
Cost
----
1.62 [m]
Free Variables
--------------
dx : 1.2 [m] Length of an element
\vec{M} : [ 1.98e+03 1.27e+03 713 317 ... ] [N*m] Internal moment
\vec{V} : [ 660 528 396 264 ... ] [N] Internal shear
\vec{\theta} : [ 0 0.177 0.285 0.341 ... ] Slope
\vec{w} : [ 0 0.106 0.384 0.759 ... ] [m] Displacement
Constants
---------
EI : 1.1e+04 [N*m**2]
L : 6 [m] Overall beam length
M_{tip} : 0 [N*m] Tip moment
V_{tip} : 0 [N] Tip loading
\theta_{base} : 0 Base angle
w_{base} : 0 [m] Base deflection
\vec{M} : [ - - - - ... ] [N*m] Internal moment
\vec{V} : [ - - - - ... ] [N] Internal shear
\vec{\theta} : [ 0 - - - ... ] Slope
\vec{q} : [ 110 110 110 110 ... ] [N/m] Distributed load at each point
\vec{w} : [ 0 - - - ... ] [m] Displacement
Sensitivities
-------------
L : 4 Overall beam length
\vec{q} : [ 0.0072 0.0416 0.118 0.234 ... ] Distributed load at each point
EI : -1
By plotting the deflection, we can see that the agreement between the analytical solution and the GP solution is good.