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: %.4g" % sol["cost"])
print("Optimal x val: %.4g" % sol["variables"][x])
Of course, the optimal value is 1. Output:
Optimal cost: 1
Optimal x val: 1
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
print(m.solve(verbosity=0).table())
The output is
┃/┓
Cost╺┫ ┃
(0.00367/m³) ┃/┣╸alpha
┃/┛ (2, fixed)
┃┓
┃┃
┃┃
┃┣╸A_{wall} = 200m²
┃┃
┃┛
┃┓
Model╺┫┃
┃┃
┃┣╸A_{wall} ≥ 2·h·w + 2·h·d
┃┃
┃┛
┃┣╸alpha = 2
┃┛
┃┣╸alpha ≤ h/w
┃┛
Free Variables
--------------
d : 8.17 [m] depth
h : 8.163 [m] height
w : 4.081 [m] width
Fixed Variables
---------------
A_{floor} : 50 [m²] upper limit, floor area
A_{wall} : 200 [m²] 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
Variable Sensitivities
----------------------
A_{wall} : -1.5 upper limit, wall area
alpha : +0.5 lower limit, wall aspect ratio
Most Sensitive Constraints
--------------------------
+1.5 : A_{wall} ≥ 2·h·w + 2·h·d
+0.5 : alpha ≤ h/w
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")
# because its units are incorrect the line below will print a warning
bad_monomial_equality = (M == V)
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.summary())
The output is:
Infeasible monomial equality: Cannot convert from 'V [m³]' to 'M [kg]'
┃┓ ┓ ┓
┃┃ ┃ ┃
┃┃ ┣╸d[0] ┣╸M
Cost╺┫┃ ┃ (0.464m) ┃ (100kg, fixed)
(1.29m²) ┃┣╸A ┛ ┛
┃┃ (1.29m²) ┓
┃┃ ┣╸d[1]·d[2]
┃┛ ┛ (0.215m²)
┃┓
┃┃
┃┣╸A ≥ 2·(d[0]·d[1] + d[0]·d[2] + d[1]·d[2])
┃┛
┃┓
┃┃
┃┣╸M = 100kg
┃┛
┃┓
Model╺┫┃
┃┣╸M = V·\rho
┃┛
┃┓
┃┃
┃┣╸V = d[0]·d[1]·d[2]
┃┛
┃┓
┃┃
┃┣╸\rho = 1,000kg/m³
┃┛
Free Variables
--------------
A : 1.293 [m²] Surface Area of the Tank
V : 0.1 [m³] Volume of the Tank
d : [ 0.464 0.464 0.464 ] [m] Dimension Vector
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 pickle
import numpy as np
from gpkit import Variable, Model, SolutionArray
pi = np.pi
# 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")
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.summary())
# save solution to a file and retrieve it
sol.save("solution.pkl")
sol.save_compressed("solution.pgz")
print(sol.diff("solution.pkl"))
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.summary())
sol_loaded = pickle.load(open("solution.pkl", "rb"))
assert sol_loaded.almost_equal(SolutionArray.decompress_file("solution.pgz"))
print(sweepsol.diff(sol_loaded, absdiff=True, senssdiff=True))
The output is:
SINGLE
======
┃┓ ┓
┃┃ ┃
┃┃ ┣╸W_0
Cost╺┫┃ ┃ (4,940N, fixed)
(303N) ┃┣╸W ┛
┃┃ (7,341N) ┓ ┓
┃┃ ┣╸W_w ┣╸S╶⎨
┃┛ ┛ (2,401N) ┛ (16.4m²)
┃┓
┃┣╸W ≥ W_0 + W_w
┃┛
┃┣╸C_D ≥ (CDA0)/S + k·C_f·(\frac{S}{S_{wet}}) + C_L²/(π·A·e)
┃┛
┃┣╸D ≥ 0.5·\rho·S·C_D·V²
┃┛
┃┣╸W_0 = 4,940N
┃┛
Model╺┫┣╸W ≤ 0.5·\rho·S·C_L·V²
┃┛
┃┣╸e = 0.95
┃┣╸(\frac{S}{S_{wet}}) = 2.05
┃┣╸C_f ≥ 0.074/Re^0.2
┃┣╸k = 1.2
┃┓
┃┃
┃┣╸[12 terms]
┃┃
┃┛
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²] total wing area
V : 38.15 [m/s] cruising speed
W : 7341 [N] total aircraft weight
W_w : 2401 [N] wing weight
Solution Diff
=============
(argument is the baseline solution)
** no constraint differences **
Relative Differences |above 1%|
-------------------------------
The largest is +0%.
SWEEP
=====
Optimal Cost
------------
[ 338 396 294 326 ]
Swept Variables
---------------
V : [ 45 55 45 55 ] [m/s] cruising speed
V_{min} : [ 20 20 25 25 ] [m/s] takeoff speed
Free Variables
--------------
A : [ 6.2 4.77 8.84 7.16 ] aspect ratio
C_D : [ 0.0146 0.0123 0.0196 0.0157 ] Drag coefficient of wing
C_L : [ 0.296 0.198 0.463 0.31 ] Lift coefficent of wing
C_f : [ 0.00333 0.00314 0.00361 0.00342 ] skin friction coefficient
D : [ 338 396 294 326 ] [N] total drag force
Re : [ 5.38e+06 7.24e+06 3.63e+06 4.75e+06 ] Reynold's number
S : [ 18.6 17.3 12.1 11.2 ] [m²] total wing area
W : [ 6.85e+03 6.4e+03 6.97e+03 6.44e+03 ] [N] total aircraft weight
W_w : [ 1.91e+03 1.46e+03 2.03e+03 1.5e+03 ] [N] wing weight
Solution Diff
=============
(argument is the baseline solution)
** no constraint differences **
Relative Differences |above 1%|
-------------------------------
Re : [ +46.4% +97.1% -1.1% +29.2% ] Reynold's number
C_L : [ -40.6% -60.2% -7.2% -37.9% ] Lift coefficent of wing
V : [ +18.0% +44.2% +18.0% +44.2% ] cruising speed
W_w : [ -20.7% -39.3% -15.6% -37.4% ] wing weight
C_D : [ -29.0% -40.4% -5.0% -23.9% ] Drag coefficient of wing
A : [ -26.7% -43.6% +4.5% -15.3% ] aspect ratio
S : [ +12.8% +5.5% -26.5% -32.0% ] total wing area
D : [ +11.5% +30.7% -2.9% +7.5% ] total drag force
V_{min} : [ -9.1% -9.1% +13.6% +13.6% ] takeoff speed
W : [ -6.8% -12.8% -5.1% -12.2% ] total aircraft weight
C_f : [ -7.3% -12.7% - -5.0% ] skin friction coefficient
Absolute Differences |above 0.1|
--------------------------------
Re : [ +1.7e+06 +3.6e+06 -4.1e+04 +1.1e+06 ] Reynold's number
W : [ -5e+02 -9.4e+02 -3.8e+02 -9e+02 ] [N] total aircraft weight
W_w : [ -5e+02 -9.4e+02 -3.8e+02 -9e+02 ] [N] wing weight
D : [ +35 +93 -8.8 +23 ] [N] total drag force
V : [ +6.8 +17 +6.8 +17 ] [m/s] cruising speed
S : [ +2.1 +0.9 -4.4 -5.3 ] [m²] total wing area
V_{min} : [ -2 -2 +3 +3 ] [m/s] takeoff speed
A : [ -2.3 -3.7 +0.38 -1.3 ] aspect ratio
C_L : [ -0.2 -0.3 - -0.19 ] Lift coefficent of wing
Sensitivity Differences |above 0.1|
-----------------------------------
V : [ +0.59 +0.97 +0.25 +0.75 ] cruising speed
V_{min} : [ -0.45 -0.67 - -0.34 ] takeoff speed
C_{L,max} : [ -0.23 -0.34 - -0.17 ] max CL with flaps down
e : [ +0.15 +0.25 - +0.19 ] Oswald efficiency factor
W_0 : [ - -0.17 - -0.16 ] aircraft weight excluding wing
\rho : [ - +0.13 - +0.19 ] density of air
(\frac{S}{S_{wet}}) : [ +0.13 +0.20 - +0.11 ] wetted area ratio
k : [ +0.13 +0.20 - +0.11 ] form factor
N_{ult} : [ -0.11 -0.18 - -0.14 ] ultimate load factor
W_{W_{coeff1}} : [ -0.11 -0.18 - -0.14 ] Wing Weight Coefficent 1
\tau : [ +0.11 +0.18 - +0.14 ] airfoil thickness to chord ratio
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 parse_variables, Model, ureg
from gpkit.small_scripts import mag
eps = 2e-4 # has to be quite large for consistent cvxopt printouts;
# normally you'd set this to something more like 1e-20
class Beam(Model):
"""Discretization of the Euler beam equations for a distributed load.
Variables
---------
EI [N*m^2] Bending stiffness
dx [m] Length of an element
L 5 [m] Overall beam length
Boundary Condition Variables
----------------------------
V_tip eps [N] Tip loading
M_tip eps [N*m] Tip moment
th_base eps [-] Base angle
w_base eps [m] Base deflection
Node Variables of length N
--------------------------
q 100*np.ones(N) [N/m] Distributed load
V [N] Internal shear
M [N*m] Internal moment
th [-] Slope
w [m] Displacement
Upper Unbounded
---------------
w_tip
"""
@parse_variables(__doc__, globals())
def setup(self, N=4):
# minimize tip displacement (the last w)
self.cost = self.w_tip = w[-1]
return {
"definition of dx": L == (N-1)*dx,
"boundary_conditions": [
V[-1] >= V_tip,
M[-1] >= M_tip,
th[0] >= th_base,
w[0] >= w_base
],
# 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 integration":
V[:-1] >= V[1:] + 0.5*dx*(q[:-1] + q[1:]),
"moment integration":
M[:-1] >= M[1:] + 0.5*dx*(V[:-1] + V[1:]),
# slope and displacement increase from base to tip (right > left)
"theta integration":
th[1:] >= th[:-1] + 0.5*dx*(M[1:] + M[:-1])/EI,
"displacement integration":
w[1:] >= w[:-1] + 0.5*dx*(th[1:] + th[:-1])
}
b = Beam(N=6, substitutions={"L": 6, "EI": 1.1e4, "q": 110*np.ones(6)})
sol = b.solve(verbosity=0)
print(sol.summary(maxcolumns=6))
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: # pragma: no cover
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:
┃┓ ┓ ┓ ┓ ┓
┃┃ ┃ ┃ ┃ ┣╸th[2]╶⎨
┃┃ ┃ ┃ ┣╸w[2] ┛ (0.285)
┃┃ ┃ ┃ ┃ (0.384m) ┣╸th[1]╶⎨
┃┃ ┃ ┣╸w[3] ┛ ┣╸w[1]╶⎨
┃┃ ┃ ┃ (0.76m) ┣╸th[3] ┣╸th[2]╶⎨
┃┃ ┃ ┃ ┛ (0.341) ┛
┃┃ ┣╸w[4] ┃ ┣╸th[2]╶⎨
┃┃ ┃ (1.18m) ┛ ┛
Cost╺┫┃ ┃ ┣╸th[3] ┣╸th[2]╶⎨
(1.62m) ┃┣╸w[5] ┃ ┛ ┛
┃┃ (1.62m) ┃ ┓ ┓
┃┃ ┃ ┣╸th[4] ┣╸th[2]╶⎨
┃┃ ┛ ┛ (0.363) ┛
┃┃ ┓ ┓
┃┃ ┣╸th[5] ┣╸th[2]╶⎨
┃┃ ┛ (0.367) ┛
┃┃ ┓ ┓
┃┃ ┣╸th[4] ┣╸th[2]╶⎨
┃┛ ┛ ┛
┃┓ ┓
┃┃ ┃
┃┃ ┣╸L = 5·dx
┃┃ ┛
┃┃ ┓
┃┃ ┃
┃┃ ┣╸L = 6m
┃┃ ┛
┃┃ ┣╸EI = 11,000N·m²
Model╺┫┃ ┣╸w[5] ≥ w[4] + 0.5·dx·(th[5] + th[4])
┃┣╸Beam ┣╸th[2] ≥ th[1] + 0.5·dx·(M[2] + M[1])/EI
┃┃ ┣╸w[4] ≥ w[3] + 0.5·dx·(th[4] + th[3])
┃┃ ┣╸M[1] ≥ M[2] + 0.5·dx·(V[1] + V[2])
┃┃ ┣╸th[3] ≥ th[2] + 0.5·dx·(M[3] + M[2])/EI
┃┃ ┣╸V[3] ≥ V[4] + 0.5·dx·(q[3] + q[4])
┃┃ ┣╸th[1] ≥ th[0] + 0.5·dx·(M[1] + M[0])/EI
┃┃ ┓
┃┃ ┃
┃┃ ┣╸[17 terms]
┃┛ ┛
Free Variables
--------------
dx : 1.2 [m] Length of an element
M : [ 1.98e+03 1.27e+03 713 317 79.2 0.0002 ] [N·m] Internal moment
V : [ 660 528 396 264 132 0.0002 ] [N] Internal shear
th : [ 0.0002 0.177 0.285 0.341 0.363 0.367 ] Slope
w : [ 0.0002 0.107 0.384 0.76 1.18 1.62 ] [m] Displacement
By plotting the deflection, we can see that the agreement between the analytical solution and the GP solution is good.