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
-------------
   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):

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.

import numpy as np
from gpkit.shortcuts import *

# Constants
k = Var("k", 1.2, "-", "form factor")
e = Var("e", 0.95, "-", "Oswald efficiency factor")
mu = Var("\\mu", 1.78e-5, "kg/m/s", "viscosity of air")
pi = Var("\\pi", np.pi, "-", "half of the circle constant")
rho = Var("\\rho", 1.23, "kg/m^3", "density of air")
tau = Var("\\tau", 0.12, "-", "airfoil thickness to chord ratio")
N_ult = Var("N_{ult}", 3.8, "-", "ultimate load factor")
V_min = Var("V_{min}", 22, "m/s", "takeoff speed")
C_Lmax = Var("C_{L,max}", 1.5, "-", "max CL with flaps down")
S_wetratio = Var("(\\frac{S}{S_{wet}})", 2.05, "-", "wetted area ratio")
W_W_coeff1 = Var("W_{W_{coeff1}}", 8.71e-5, "1/m", "Wing Weight Coefficent 1")
W_W_coeff2 = Var("W_{W_{coeff2}}", 45.24, "Pa", "Wing Weight Coefficent 2")
CDA0 = Var("(CDA0)", 0.031, "m^2", "fuselage drag area")
W_0 = Var("W_0", 4940.0, "N", "aircraft weight excluding wing")

# Free Variables
D = Var("D", "N", "total drag force")
A = Var("A", "-", "aspect ratio")
S = Var("S", "m^2", "total wing area")
V = Var("V", "m/s", "cruising speed")
W = Var("W", "N", "total aircraft weight")
Re = Var("Re", "-", "Reynold's number")
C_D = Var("C_D", "-", "Drag coefficient of wing")
C_L = Var("C_L", "-", "Lift coefficent of wing")
C_f = Var("C_f", "-", "skin friction coefficient")
W_w = Var("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=1)

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)
m.solve(verbosity=1)

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
import gpkit
from gpkit.shortcuts import Var, Vec, Model
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 setup(self, N=4):
        EI = Var("EI", 1e4, "N*m^2")
        dx = Var("dx", "m", "Length of an element")
        L = Var("L", 5, "m", "Overall beam length")
        q = Vec(N, "q", 100*np.ones(N), "N/m", "Distributed load at each point")
        V = Vec(N, "V", "N", "Internal shear")
        V_tip = Var("V_{tip}", 0, "N", "Tip loading")
        M = Vec(N, "M", "N*m", "Internal moment")
        M_tip = Var("M_{tip}", 0, "N*m", "Tip moment")
        th = Vec(N, "\\theta", "-", "Slope")
        th_base = Var("\\theta_{base}", 0, "-", "Base angle")
        w = Vec(N, "w", "m", "Displacement")
        w_base = Var("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)
        return w[-1], [shear_eq, moment_eq, theta_eq, displ_eq, L == (N-1)*dx]


b = Beam(N=10, substitutions={"L": 6, "EI": 1.1e4, "q": 110*np.ones(10)})
sol = b.solve(verbosity=1)
w_gp = sol("w")  # deflection along beam

L, EI, q = sol("L"), sol("EI"), sol("q")
x = np.linspace(0, mag(L), len(q))*gpkit.units.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)) <= 1e-2*gpkit.units.m

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
----
 1.473 [m] 

Free Variables
--------------
    dx : 0.6667                                          [m]   Length of an element
     M : [ 1.8e+03   1.42e+03  1.09e+03  800      ... ]  [N*m] Internal moment     
     V : [ 600       533       467       400      ... ]  [N]   Internal shear      
\theta : [  -        0.0976    0.174     0.231    ... ]        Slope               
     w : [  -        0.0325    0.123     0.258    ... ]  [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               
            M : [  -         -         -         -       ... ]  [N*m]    Internal moment               
            V : [  -         -         -         -       ... ]  [N]      Internal shear                
       \theta : [ 0          -         -         -       ... ]           Slope                         
            q : [ 100       100       100       100      ... ]  [N/m]    Distributed load at each point
            w : [ 0          -         -         -       ... ]  [m]      Displacement                  

Sensitivities
-------------
 L : 4                                              Overall beam length           
 q : [ 0.0013    0.00762   0.0223    0.0454   ... ] 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.