Examples¶
iPython Notebook Examples¶
More examples, including some in-depth and experimental models, 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\).
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.
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=1)
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
-------------
A_{wall} : -1.5 upper limit, wall area
alpha : 0.5 lower limit, wall aspect ratio
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):
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=1)
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
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.
The output is
Using solver 'cvxopt'
Solving for 9 variables.
Solving took 0.0492 seconds.
Cost
----
255
Free Variables
--------------
A : 12.7 Aspect ratio
C_D : 0.0231 Drag coefficient
C_L : 0.6512 Lift coefficient
C_f : 0.003857 Skin friction coefficient
Re : 2.598e+06 Reynolds number
S : 12.08 [m**2] Wing planform area
V : 38.55 [m/s] Cruise velocity
W : 7189 [N] Total aircraft weight
W_w : 2249 [N] Wing weight
Constants
---------
(CDA_0) : 0.0306 [m**2] Fuselage drag area
C_{L,max} : 2 Maximum C_L, flaps down
N_{lift} : 2.5 Ultimate load factor
V_{min} : 22 [m/s] Desired landing speed
W_0 : 4940 [N] Aircraft weight excluding wing
\frac{S_{wet}}{S} : 2.05 Wetted area ratio
\mu : 1.78e-05 [kg/m/s] Viscosity of air
\rho : 1.23 [kg/m**3] Density of air
\tau : 0.12 Airfoil thickness-to-chord ratio
cww1 : 45.42 [N/m**2] Wing weight area factor
cww2 : 8.71e-05 [1/m] Wing weight bending factor
e : 0.96 Oswald efficiency factor
k : 1.2 Form factor
Sensitivities
-------------
(CDA_0) : 0.1097 Fuselage drag area
C_{L,max} : -0.1307 Maximum C_L, flaps down
N_{lift} : 0.2923 Ultimate load factor
V_{min} : -0.2614 Desired landing speed
W_0 : 0.9953 Aircraft weight excluding wing
\frac{S_{wet}}{S} : 0.4108 Wetted area ratio
\mu : 0.08217 Viscosity of air
\rho : -0.1718 Density of air
\tau : -0.2923 Airfoil thickness-to-chord ratio
cww1 : 0.09428 Wing weight area factor
cww2 : 0.2923 Wing weight bending factor
e : -0.4795 Oswald efficiency factor
k : 0.4108 Form factor
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.shortcuts import *
def beam(N=10, L=5., EI=1E4, P=100):
dx = Var("dx", L/(N-1), units="m")
EI = Var("EI", EI, units="N*m^2")
p = Vec(N, "p", units="N/m", label="Distributed load")
p = p.sub(p, P*np.ones(N))
V = Vec(N, "V", units="N", label="Internal shear")
M = Vec(N, "M", units="N*m", label="Internal moment")
th = Vec(N, "th", units="-", label="Slope")
w = Vec(N, "w", units="m", label="Displacement")
eps = 1E-16 #an arbitrarily small positive number
substitutions = {var: eps for var in [V[-1], M[-1], th[0], w[0]]}
objective = w[-1]
constraints = [V.left[1:N] >= V[1:N] + 0.5*dx*(p.left[1:N] + p[1:N]),
M.left[1:N] >= M[1:N] + 0.5*dx*(V.left[1:N] + V[1:N]),
th.right[0:N-1] >= th[0:N-1] + 0.5*dx*(M.right[0:N-1] + M[0:N-1])/EI,
w.right[0:N-1] >= w[0:N-1] + 0.5*dx*(th.right[0:N-1]+ th[0:N-1])
]
return Model(objective, constraints, substitutions)
N = 10 # [-] grid size
L = 5. # [m] beam length
EI = 1E4 # [N*m^2] elastic modulus * area moment of inertia
P = 100 # [N/m] magnitude of distributed load
m = beam(N, L, EI, P)
sol = m.solve(verbosity=1)
x = np.linspace(0, L, N) # position along beam
w_gp = sol("w") # deflection along beam
w_exact = P/(24.*EI)* x**2 * (x**2 - 4*L*x + 6*L**2) # analytical soln
assert max(abs(w_gp - w_exact)) <= 1e-2
PLOT = False
if PLOT:
import matplotlib.pyplot as plt
x_exact = np.linspace(0, L, 1000)
w_exact = P/(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
----
0.7813 [m]
Free Variables
--------------
M : [ 1.25e+03 988 756 556 ... ] [N*m] Internal moment
V : [ 500 444 389 333 ... ] [N] Internal shear
th : [ - 0.0622 0.111 0.147 ... ] Slope
w : [ - 0.0173 0.0653 0.137 ... ] [m] Displacement
Constants
---------
EI : 1e+04 [N*m**2]
dx : 0.5556 [m]
M : [ - - - - ... ] [N*m] Internal moment
V : [ - - - - ... ] [N] Internal shear
th : [ 1e-16 - - - ... ] Slope
w : [ 1e-16 - - - ... ] [m] Displacement
Sensitivities
-------------
EI : -1
dx : 4
By plotting the deflection, we can see that the agreement between the analytical solution and the GP solution is good.
