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.

_images/beam.svg