GPkit Documentation¶

Table of contents¶

Geometric Programming 101¶

What is a GP?¶

A Geometric Program (GP) is a special type of constrained non-linear optimization problem.

A GP is made up of special types of functions called monomials and posynomials. In the context of a GP, a monomial is defined as:

\[f(x) = c x_1^{a_1} x_2^{a_2} ... x_n^{a_n}\]

where \(c\) is a positive constant, \(x_{1..n}\) are the decision variables, and \(a_{1..n}\) are real exponents.

Building on this, a posynomial is defined as a sum of monomials:

\[g(x) = \sum_{k=1}^K c_k x_1^{a_1k} x_2^{a_2k} ... x_n^{a_nk}\]

Using these definitions, a GP in Standard Form is written as:

\[\begin{split}\begin{array}[lll]\text{} \text{minimize} & f_0(x) & \\ \text{subject to} & f_i(x) = 1, & i = 1,....,m \\ & g_i(x) \leq 1, & i = 1,....,n \end{array}\end{split}\]

Why are GPs special?¶

Geometric programs have several powerful properties:

1. Unlike most non-linear optimization problems, large GPs can be solved extremely quickly.
2. If there exists an optimal solution to a GP, it is guaranteed to be globally optimal.
3. Modern GP solvers require no initial guesses or tuning of solver parameters.

These properties arise because GPs become convex optimization problems via a logarithmic transformation. In addition to their mathematical benefits, recent research has shown that many practical problems can be formulated as GPs or closely approximated as GPs.

Where can I learn more?¶

To learn more about GPs, refer to the following resources:

• A tutorial on geometric programming, by S. Boyd, S.J. Kim, L. Vandenberghe, and A. Hassibi.
• Convex optimization, by S. Boyd and L. Vandenberghe.

Installation Instructions¶

If you encounter any bugs during installation, please email gpkit@mit.edu.

Mac OS X¶

1. Install Python and build dependencies¶

• Install the Python 2.7 version of Anaconda.
• If you don’t want to install Anaconda, you’ll need gcc and pip, and may find sympy, scipy, and iPython Notebook useful.
• If which gcc does not return anything, install the Apple Command Line Tools.
• _Optional:_ to install gpkit into an isolated python environment you can create a new conda virtual environment with conda create -n gpkit anaconda and activate it with source activate gpkit.
• Run pip install ctypesgen in the Terminal.

2. Install either the MOSEK or CVXOPT GP solvers¶

• Download CVXOPT, then:

  • Read the official instructions and requirements
  • In the Terminal, navigate to the cvxopt folder
  • Run python setup.py install

• Download MOSEK, then:

  • Move the mosek folder to your home directory
  • Follow these steps for Mac.
  • Request an academic license file and put it in ~/mosek/

3. Install GPkit¶

• Run pip install gpkit at the command line.
• Run python -c "import gpkit.tests; gpkit.tests.run()"
• If you want units support, install pint with pip install pint.

Linux¶

1. Install either the MOSEK or CVXOPT GP solvers¶

• Download CVXOPT, then:

  • Read the official instructions and requirements
  • In a terminal, navigate to the cvxopt folder
  • Run python setup.py install

• Download MOSEK, then:

  • Move the mosek folder to your home directory
  • Follow these steps for Linux.
  • Request an academic license file and put it in ~/mosek/

2. Install GPkit¶

• _Optional:_ to install gpkit into an isolated python environment, install virtualenv, run virtualenv $DESTINATION_DIR then activate it with source activate $DESTINATION_DIR/bin.
• Run pip install ctypesgen at the command line.
• Run pip install gpkit at the command line.
• Run python -c "import gpkit.tests; gpkit.tests.run()"
• If you want units support, install pint with pip install pint.
• You may find sympy, scipy, and iPython Notebook to be useful additional packages as well.

Windows¶

1. Install Python dependencies¶

• Install the Python 2.7 version of Anaconda.
• If you don’t want to install Anaconda, you’ll need gcc and pip, and may find sympy, scipy, and iPython Notebook useful.
• _Optional:_ to install gpkit into an isolated python environment you can create a new conda virtual environment with conda create -n gpkit anaconda and activate it with source activate gpkit.
• Run pip install ctypesgen at an Anaconda Command Prompt.

2. Install either the MOSEK or CVXOPT GP solvers¶

• Download CVXOPT, then follow these steps to install a linear algebra library

• Download MOSEK, then:

  • Follow these steps for Windows.
  • Request an academic license file and put it in ~/mosek/

3. Install GPkit¶

• Run pip install gpkit at an Anaconda Command Prompt.
• Run python -c "import gpkit.tests; gpkit.tests.run()"
• If you want units support, install pint with pip install pint.

Getting Started¶

GPkit is a Python package. We assume basic familiarity with Python. If you are new to Python take a look at Learn Python.

GPkit is also a command line tool. This means that you need to be in the terminal (OS X/Linux) or command prompt (Windows) to use it. If you are not familiar with working in the command line, check out this Learn Code the Hard Way tutorial.

The first thing to do is install GPkit . Once you have done this, you can start using GPkit in 3 easy steps:

1. Open your command line interface (terminal/Command Prompt).
2. Open a Python interpreter. This can be done by typing python (or ipython if installed).
3. Type import gpkit.

After doing this, your command line will look something like the following:

$ python
>>> import gpkit
>>>

From here, you can use GPkit commands to formulate and solve geometric programs. To learn how, see Basic Commands.

Writing GPkit Scripts¶

Another way to use GPkit is to write a script and save it as a .py file. To run this file (e.g. myscript.py), type the following in your command line:

$ python myscript.py

Again, ipython will also work here.

Basic Commands¶

Importing Modules¶

The first thing to do when using GPkit is to import the classes and modules you will need. For example,

from gpkit import Variable, VectorVariable, GP

Declaring Variables¶

Instances of the Variable class represent scalar decision variables. They store a key (i.e. name) used to look up the Variable in dictionaries, and optionally units, a description, and a value (if the Variable is to be held constant).

Decision Variables¶

# Declare a variable, x
x = Variable('x')

# Declare a variable, y, with units of meters
y = Variable('y','m')

# Declare a variable, z, with units of meters, and a description
z = Variable('z', 'm', 'A variable called z with units of meters')

Note: make sure you have imported the class Variable beforehand.

Fixed Variables¶

To declare a variable with a constant value, use the Variable class, as above, but specify the value= input argument:

# Declare '\\rho' equal to 1.225 kg/m^3
rho = Variable('\\rho', 1.225, 'kg/m^3', 'Density of air at sea level')

In the example above, the key name '\\rho' is for LaTeX printing (described later). The unit and description arguments are optional.

#Declare pi equal to 3.14
pi = Variable('\\pi', 3.14)

Vector Variables¶

Vector variables are represented by the VectorVariable class. The first argument is the length of the vector. All other inputs follow those of the Variable class.

# Declare a 3-element vector variable 'x' with units of 'm'
x = VectorVariable(3, "x", "m", "3-D Position")

Creating Monomials and Posynomials¶

Monomial and posynomial expressions can be created using mathematical operations on variables. This is implemented under-the-hood using operator overloading in Python.

# create a Monomial term xy^2/z
x = Variable('x')
y = Variable('y')
z = Variable('z')
m = x * y**2 / z
type(m)  # gpkit.nomials.Monomial

# create a Posynomial expression x + xy^2
x = Variable('x')
y = Variable('y')
p = x + x * y**2
type(p)  # gpkit.nomials.Posynomial

Declaring Constraints¶

Constraint objects represent constraints of the form Monomial >= Posynomial or Monomial == Monomial (which are the forms required for GP-compatibility).

Note that constraints must be formed using <=, >=, or == operators, not < or >.

# consider a block with dimensions x, y, z less than 1
# constrain surface area less than 1.0 m^2
x = Variable('x', 'm')
y = Variable('y', 'm')
z = Variable('z', 'm')
S = Variable('S', 1.0, 'm^2')
c = (2*x*y + 2*x*z + 2*y*z <= S)
type(c)  # gpkit.nomials.Constraint

Declaring Objective Functions¶

To declare an objective function, assign a Posynomial (or Monomial) to a variable name, such as objective.

objective = 1/(x*y*z)

By convention, the objective is the function to be minimized. If you wish to maximize a function, take its reciprocal. For example, the code above creates an objective which, when minimized, will maximize x*y*z.

Formulating a GP¶

The GP class represents a geometric programming problem. To create one, pass an objective and list of Constraints:

objective = 1/(x*y*z)
constraints = [2*x*y + 2*x*z + 2*y*z <= S,
               x >= 2*y]
gp = GP(objective, constraints)

Solving the GP¶

sol = gp.solve()

Printing Results¶

print sol.table()

print "The x dimension is %s." % (sol(x))

Examples¶

A Trivial GP¶

from gpkit import Variable, GP

# Decision variable
x = Variable('x')

# Constraint
constraints = [x >= 1]

# Objective (to minimize)
objective = x

# Formulate the GP
gp = GP(objective, constraints)

# Solve the GP
sol = gp.solve()

# Print results table
print sol.table()

Maximizing the Volume of a Box¶

from gpkit import Variable, GP

# 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 GP
gp = GP(objective, constraints)

# Solve the GP
sol = gp.solve()

# Print results table
print sol.table()

Water Tank¶

from gpkit import Variable, VectorVariable, GP
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
               )

gp = GP(A, constraints)
sol = gp.solve(printing=False)
print sol(A)
print sol(V)
print sol(d)

iPython Notebook Examples¶

Also available on nbviewer.

BOX¶

Maximize volume while limiting the aspect ratios and area of the sides.

This is a simple demonstration of gpkit, based on an example from A Tutorial on Geometric Programming by Boyd et al..

Set up the modelling environment¶

First we’ll to import GPkit and turn on \(\LaTeX\) printing for GPkit variables and equations.

import gpkit
import gpkit.interactive
gpkit.interactive.init_printing()

Now we declare the optimisation parameters:

alpha = gpkit.Variable("\\alpha", 2, "-", "lower limit, wall aspect ratio")
beta = gpkit.Variable("\\beta", 10, "-", "upper limit, wall aspect ratio")
gamma = gpkit.Variable("\\gamma", 2, "-", "lower limit, floor aspect ratio")
delta = gpkit.Variable("\\delta", 10, "-", "upper limit, floor aspect ratio")
A_wall = gpkit.Variable("A_{wall}", 200, "m^2", "upper limit, wall area")
A_floor = gpkit.Variable("A_{floor}", 50, "m^2", "upper limit, floor area")

Next, we declare the decision variables:

var = gpkit.Variable # a convenient shorthand

h = var("h", "m", "height")
w = var("w", "m", "width")
d = var("d", "m", "depth")

Then we form equations of the system:

V = h*w*d
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]

Formulate the optimisation problem¶

Here we write the optimisation problem as a standard form geometric program. Note that by putting \(\tfrac{1}{V}\) as our cost function, we are maximizing \(V\).

gp = gpkit.GP(1/V, constraints)

Now we can check that our equations are correct by using the built-in latex printing.

gp

\[\begin{split}\begin{array}[ll] \text{} \text{minimize} & \frac{1}{d h w}\mathrm{\left[ \tfrac{1}{m^{3}} \right]} \\ \text{subject to} & A_{wall} \geq 2d h + 2h w \\ & A_{floor} \geq d w \\ & \alpha \leq \frac{h}{w} \\ & \beta \geq \frac{h}{w} \\ & \gamma \leq \frac{d}{w} \\ & \delta \geq \frac{d}{w} \\ \text{substituting} & A_{floor} = 50 \\ & A_{wall} = 200 \\ & \alpha = 2 \\ & \beta = 10 \\ & \delta = 10 \\ & \gamma = 2 \\ \end{array}\end{split}\]

That looks fine, but let’s change \(A_{floor}\) to \(100\), just for fun.

Note that replace=True will be needed, since \(A_{floor}\) has already been substituted in.

gp.sub(A_floor, 100, replace=True) # (var, val) substitution syntax

And heck, why not change \(\alpha\) and \(\gamma\) to \(1\) while we’re at it?

gp.sub({alpha: 1, gamma: 1}, replace=True) # ({var1: val1, var2: val2}) substitution syntax

Now check that those changes took:

gp

\[\begin{split}\begin{array}[ll] \text{} \text{minimize} & \frac{1}{d h w}\mathrm{\left[ \tfrac{1}{m^{3}} \right]} \\ \text{subject to} & A_{wall} \geq 2d h + 2h w \\ & A_{floor} \geq d w \\ & \alpha \leq \frac{h}{w} \\ & \beta \geq \frac{h}{w} \\ & \gamma \leq \frac{d}{w} \\ & \delta \geq \frac{d}{w} \\ \text{substituting} & A_{floor} = 100 \\ & A_{wall} = 200 \\ & \alpha = 1 \\ & \beta = 10 \\ & \delta = 10 \\ & \gamma = 1 \\ \end{array}\end{split}\]

Looks good!

Solve the GP¶

sol = gp.solve()

Using solver 'cvxopt'
Solving took 0.0184 seconds

Analyse results¶

print sol.table()

0.0025981 : Cost
          | Free variables
        d : 11.5     [m] depth
        h : 5.77     [m] height
        w : 5.77     [m] width
          |
          | Constants
A_{floor} : 100      [m**2] upper limit, floor area
 A_{wall} : 200      [m**2] upper limit, wall area
   alpha : 1        [-] lower limit, wall aspect ratio
    beta : 10       [-] upper limit, wall aspect ratio
   delta : 10       [-] upper limit, floor aspect ratio
   gamma : 1        [-] lower limit, floor aspect ratio
          |
          | Constant sensitivities
 A_{wall} : -1.5     [-] upper limit, wall area
   alpha : 0.5      [-] lower limit, wall aspect ratio
          |

Hmm, why didn’t \(A_{floor}\) show up in the sensitivities list?

sol.senssubinto(A_floor)

-3.4698494246624108e-09

Its sensitivity is tiny; changing it near this value doesn’t affect the cost at all; that constraint is loose! Let’s sweep over a range of \(A_{floor}\) values to figure out where it becomes loose.

Sweep and plot results¶

Import the plotting library matplotlib and the math library numpy:

%matplotlib inline
%config InlineBackend.figure_format = 'retina' # for high-DPI displays
import matplotlib.pyplot as plt
import numpy as np

Solve for values of \(A_{floor}\) from 10 to 100 using a “sweep” substitution:

gp.sub(A_floor, ('sweep', np.linspace(10, 100, 50)), replace=True)
sol = gp.solve()
print sol.table()

Using solver 'cvxopt'
Sweeping 1 variables over 50 passes
Solving took 0.547 seconds
 0.0032211 : Cost (average of 50 values)
           | Free variables (average)
         d : 8.33     [m] depth
         h : 7.74     [m] height
         w : 5.69     [m] width
           |
           | Constants (average)
 A_{floor} : 55       [m**2] upper limit, floor area
  A_{wall} : 200      [m**2] upper limit, wall area
    alpha : 1        [-] lower limit, wall aspect ratio
     beta : 10       [-] upper limit, wall aspect ratio
    delta : 10       [-] upper limit, floor aspect ratio
    gamma : 1        [-] lower limit, floor aspect ratio
           |
           | Constant sensitivities (average)
 A_{floor} : -0.274   [-] upper limit, floor area
  A_{wall} : -1.23    [-] upper limit, wall area
    alpha : 0.226    [-] lower limit, wall aspect ratio
           |

It seems we got some sensitivity out of \(A_{floor}\) on average over these points; let’s plot it:

plt.plot(sol(A_floor), sol(d), linewidth=1, alpha=0.5)
plt.plot(sol(A_floor), sol(h), linewidth=1, alpha=0.5)
plt.plot(sol(A_floor), sol(w), '--', linewidth=2, alpha=0.5)
plt.legend(['depth', 'height', 'width'])
plt.ylabel("Optimal dimensions [m]")
_ = plt.xlabel("A_floor [m^2]") # the _ catches the returned label object, since we don't need it

There’s an interesting elbow when \(A_{floor} \approx 50 \mathrm{\ m^2}\).

Interactive analysis¶

Let’s investigate it with the cadtoons library. Running cadtoon.py box.svg in this folder creates an interactive SVG graphic for us.

First, import the functions to display HTML in iPython Notebook, and the ractivejs library.

from IPython.display import HTML, display
from string import Template

Then write controls that link the optimization to the animation. It’s generally very helpful to play around with with writing constraints in the cadtoons sandbox before copying the javascript code from there to iPython.

ractor = Template("""
var w = $w/10,
    h = $h/10,
    d = $d/10""")

def ractorfn(sol):
    return ractor.substitute(w=sol(w), d=sol(d), h=sol(h))

constraints="""
var dw = 50*(w-1),
    dh = 50*(h-1),
    dd = 50*(d-1)

r.plan.scalex = w
r.plan.scaley = d

r.top.scalex = w
r.top.scaley = h
r.top.y = -dh -dd
r.bottom.scalex = w
r.bottom.scaley = h
r.bottom.y = dh + dd

r.left.scalex = h
r.left.scaley = d
r.left.x = -dh - dw
r.right.scalex = h
r.right.scaley = d
r.right.x = dh + dw """

def ractivefn(gp):
    sol = gp.solution
    live = "<script>" + ractorfn(sol) + constraints + "</script>"
    display(HTML(live))
    # if you enable the line below, you can try navigating the sliders by sensitivities
    # print sol.table(["cost", "sensitivities"])

Now that the display window has been created, let’s use an iPython widget to explore the “design space” of our box. This widget runs the gpkit solver every time a slider is changed, so you only solve the points you see.

with open("box.gpkit", 'r') as file:
    display(HTML(file.read()))
display(HTML(display(HTML("<style> #ractivecontainer"
             "{position:absolute; height: 0;"
             "right: 0; top: 5em;} </style>"))))

<IPython.core.display.HTML at 0x7fd2dce97850>

gpkit.interactive.widget(gp, ractivefn, {
    "A_{floor}": (10, 1000, 10), "A_{wall}": (10, 1000, 10),
    "\\delta": (0.1, 20, 0.1), "\\gamma": (0.1, 20, 0.1),
    "\\alpha": (0.1, 20, 0.1), "\\beta": (0.1, 20, 0.1), })

Now that we’ve found a good range to explore, we’ll make and interactive javascript widget that presolves all options and works in a iPython notebook hosted by nbviewer.

gpkit.interactive.jswidget(gp, ractorfn, constraints, {
    "A_{floor}": (10, 100, 10), "A_{wall}": (10, 210, 20),
    "\\delta": (1, 10, 3), "\\gamma": (1, 10, 3),
    "\\alpha": (1, 10, 3), "\\beta": (1, 10, 3), })

This concludes the Box example. Try playing around with the sliders up above until you’re bored of this simple system; then check out one of the other examples. Thanks for reading!

Import CSS for nbviewer¶

If you have a local iPython stylesheet installed, this will add it to the iPython Notebook:

from IPython import utils
from IPython.core.display import HTML
import os
def css_styling():
    """Load default custom.css file from ipython profile"""
    base = utils.path.get_ipython_dir()
    styles = ("<style>\n%s\n</style>" %
              open(os.path.join(base,'profile_default/static/custom/custom.css'),'r').read())
    return HTML(styles)
css_styling()

AIRPLANE FUEL¶

Minimize fuel needed for a plane that can sprint and land quickly.

Set up the modelling environment¶

First we’ll to import GPkit and turn on \(\LaTeX\) printing for GPkit variables and equations.

import numpy as np
import gpkit
import gpkit.interactive
gpkit.interactive.init_printing()

declare constants¶

mon = gpkit.Variable
vec = gpkit.VectorVariable

N_lift         = mon("N_{lift}", 6.0, "-", "Wing loading multiplier")
pi             = mon("\\pi", np.pi, "-", "Half of the circle constant")
sigma_max      = mon("\\sigma_{max}", 250e6, "Pa", "Allowable stress, 6061-T6")
sigma_maxshear = mon("\\sigma_{max,shear}", 167e6, "Pa", "Allowable shear stress")
g              = mon("g", 9.8, "m/s^2", "Gravitational constant")
w              = mon("w", 0.5, "-", "Wing-box width/chord")
r_h            = mon("r_h", 0.75, "-", "Wing strut taper parameter")
f_wadd         = mon("f_{wadd}", 2, "-", "Wing added weight fraction")
W_fixed        = mon("W_{fixed}", 14.7e3, "N", "Fixed weight")
C_Lmax         = mon("C_{L,max}", 1.5, "-", "Maximum C_L, flaps down")
rho            = mon("\\rho", 0.91, "kg/m^3", "Air density, 3000m")
rho_sl         = mon("\\rho_{sl}", 1.23, "kg/m^3", "Air density, sea level")
rho_alum       = mon("\\rho_{alum}", 2700, "kg/m^3", "Density of aluminum")
mu             = mon("\\mu", 1.69e-5, "kg/m/s", "Dynamic viscosity, 3000m")
e              = mon("e", 0.95, "-", "Wing spanwise efficiency")
A_prop         = mon("A_{prop}", 0.785, "m^2", "Propeller disk area")
eta_eng        = mon("\\eta_{eng}", 0.35, "-", "Engine efficiency")
eta_v          = mon("\\eta_v", 0.85, "-", "Propeller viscous efficiency")
h_fuel         = mon("h_{fuel}", 42e6, "J/kg", "fuel heating value")
V_sprint_reqt  = mon("V_{sprintreqt}", 150, "m/s", "sprint speed requirement")
W_pay          = mon("W_{pay}", 500*9.81, "N")
R_min          = mon("R_{min}", 1e6, "m", "Minimum airplane range")
V_stallmax     = mon("V_{stall,max}", 40, "m/s", "Stall speed")
# sweep variables
R_min          = mon("R_{min}", np.linspace(1e6, 1e7, 10), "m", "Minimum airplane range")
V_stallmax     = mon("V_{stall,max}", np.linspace(30, 50, 10), "m/s", "Stall speed")

declare free variables¶

V        = vec(3, "V", "m/s", "Flight speed")
C_L      = vec(3, "C_L", "-", "Wing lift coefficent")
C_D      = vec(3, "C_D", "-", "Wing drag coefficent")
C_Dfuse  = vec(3, "C_{D_{fuse}}", "-", "Fuselage drag coefficent")
C_Dp     = vec(3, "C_{D_p}", "-", "drag model parameter")
C_Di     = vec(3, "C_{D_i}", "-", "drag model parameter")
T        = vec(3, "T", "N", "Thrust force")
Re       = vec(3, "Re", "-", "Reynold's number")
W        = vec(3, "W", "N", "Aircraft weight")
eta_i    = vec(3, "\\eta_i", "-", "Aircraft efficiency")
eta_prop = vec(3, "\\eta_{prop}", "-")
eta_0    = vec(3, "\\eta_0", "-")
W_fuel   = vec(2, "W_{fuel}", "N", "Fuel weight")
z_bre    = vec(2, "z_{bre}", "-")
S        = mon("S", "m^2", "Wing area")
R        = mon("R", "m", "Airplane range")
A        = mon("A", "-", "Aspect Ratio")
I_cap    = mon("I_{cap}", "m^4", "Spar cap area moment of inertia per unit chord")
M_rbar   = mon("\\bar{M}_r", "-")
P_max    = mon("P_{max}", "W")
V_stall  = mon("V_{stall}", "m/s")
nu       = mon("\\nu", "-")
p        = mon("p", "-")
q        = mon("q", "-")
tau      = mon("\\tau", "-")
t_cap    = mon("t_{cap}", "-")
t_web    = mon("t_{web}", "-")
W_cap    = mon("W_{cap}", "N")
W_zfw    = mon("W_{zfw}", "N", "Zero fuel weight")
W_eng    = mon("W_{eng}", "N")
W_mto    = mon("W_{mto}", "N", "Maximum takeoff weight")
W_pay    = mon("W_{pay}", "N")
W_tw     = mon("W_{tw}", "N")
W_web    = mon("W_{web}", "N")
W_wing   = mon("W_{wing}", "N")

Let’s check that the vector constraints are working:

W == 0.5*rho*C_L*S*V**2

\[\begin{split}\begin{bmatrix}{W}_{0} = 0.5S \rho {V}_{0}^{2} {C_L}_{0} & {W}_{1} = 0.5S \rho {V}_{1}^{2} {C_L}_{1} & {W}_{2} = 0.5S \rho {V}_{2}^{2} {C_L}_{2}\end{bmatrix}\end{split}\]

Looks good!

Form the optimization problem¶

In the 3-element vector variables, indexs 0, 1, and 2 are the outbound, return and sprint flights.

steady_level_flight = (W == 0.5*rho*C_L*S*V**2,
                       T >= 0.5*rho*C_D*S*V**2,
                       Re == (rho/mu)*V*(S/A)**0.5)

landing_fc = (W_mto <= 0.5*rho_sl*V_stall**2*C_Lmax*S,
              V_stall <= V_stallmax)

sprint_fc = (P_max >= T[2]*V[2]/eta_0[2],
             V[2] >= V_sprint_reqt)

drag_model = (C_D >= (0.05/S)*gpkit.units.m**2 +C_Dp + C_L**2/(pi*e*A),
              1 >= (2.56*C_L**5.88/(Re**1.54*tau**3.32*C_Dp**2.62) +
                   3.8e-9*tau**6.23/(C_L**0.92*Re**1.38*C_Dp**9.57) +
                   2.2e-3*Re**0.14*tau**0.033/(C_L**0.01*C_Dp**0.73) +
                   1.19e4*C_L**9.78*tau**1.76/(Re*C_Dp**0.91) +
                   6.14e-6*C_L**6.53/(Re**0.99*tau**0.52*C_Dp**5.19)))

propulsive_efficiency = (eta_0 <= eta_eng*eta_prop,
                         eta_prop <= eta_i*eta_v,
                         4*eta_i + T*eta_i**2/(0.5*rho*V**2*A_prop) <= 4)

# 4th order taylor approximation for e^x
z_bre_sum = 0
for i in range(1,5):
    z_bre_sum += z_bre**i/np.math.factorial(i)

range_constraints = (R >= R_min,
                     z_bre >= g*R*T[:2]/(h_fuel*eta_0[:2]*W[:2]),
                     W_fuel/W[:2] >= z_bre_sum)

weight_relations = (W_pay >= 500*g*gpkit.units.kg,
                    W_tw >= W_fixed + W_pay + W_eng,
                    W_zfw >= W_tw + W_wing,
                    W_eng >= 0.0372*P_max**0.8083 * gpkit.units.parse_expression('N/W^0.8083'),
                    W_wing/f_wadd >= W_cap + W_web,
                    W[0] >= W_zfw + W_fuel[1],
                    W[1] >= W_zfw,
                    W_mto >= W[0] + W_fuel[0],
                    W[2] == W[0])

wing_structural_model = (2*q >= 1 + p,
                         p >= 1.9,
                         tau <= 0.15,
                         M_rbar >= W_tw*A*p/(24*gpkit.units.N),
                         .92**2/2.*w*tau**2*t_cap >= I_cap * gpkit.units.m**-4 + .92*w*tau*t_cap**2,
                         8 >= N_lift*M_rbar*A*q**2*tau/S/I_cap/sigma_max * gpkit.units.parse_expression('Pa*m**6'),
                         12 >= A*W_tw*N_lift*q**2/tau/S/t_web/sigma_maxshear,
                         nu**3.94 >= .86*p**-2.38 + .14*p**0.56,
                         W_cap >= 8*rho_alum*g*w*t_cap*S**1.5*nu/3/A**.5,
                         W_web >= 8*rho_alum*g*r_h*tau*t_web*S**1.5*nu/3/A**.5
                         )

eqns = (weight_relations + range_constraints + propulsive_efficiency
        + drag_model + steady_level_flight + landing_fc + sprint_fc + wing_structural_model)

gp = gpkit.GP(W_fuel.sum(), eqns)

Design a hundred airplanes¶

sol = gp.solve()

Using solver 'cvxopt'
Sweeping 2 variables over 100 passes
Solving took 9.96 seconds

The “local model” is the power-law tangent to the Pareto frontier, gleaned from sensitivities.

sol["local_model"][0].sub(gp.substitutions)

\[0.006484\frac{R_{min}}{V_{stall,max}^{0.59}}\]

plot design frontiers¶

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

plot_frontiers = gpkit.interactive.plot_frontiers
plot_frontiers(gp, [V[0], V[1], V[2]])
plot_frontiers(gp, [S, W_zfw, P_max])
plot_frontiers(gp, ['S{\\rho_{alum}}', 'S{h_{fuel}}', 'S{A_{prop}}'])

Interactive analysis¶

Let’s investigate it with the cadtoons library. Running cadtoon.py flightconditions.svg in this folder creates an interactive SVG graphic for us.

First, import the functions to display HTML in iPython Notebook, and the ractivejs library.

from IPython.display import HTML, display
from string import Template

ractor = Template("""
var W_eng = $W_eng,
    lam = $lam

r.shearinner.scalex = 1-$tcap*10
r.shearinner.scaley = 1-$tweb*100
r.airfoil.scaley = $tau/0.13
r.fuse.scalex = $W_fus/24000
r.wing.scalex = $b/2/14
r.wing.scaley = $cr*1.21
""")

def ractorfn(sol):
    return ractor.substitute(lam = 0.5*(sol(p) - 1),
                             b = sol((S*A)**0.5),
                             cr = sol(2/(1+0.5*(sol(p) - 1))*(S/A)**0.5),
                             tcap = sol(t_cap/tau),
                             tweb = sol(t_web/w),
                             tau = sol(tau),
                             W_eng = sol(W_eng),
                             W_fus = sol(W_mto) - sol(W_wing) - sol(W_eng))

constraints = """
r.engine1.scale = Math.pow(W_eng/3000, 2/3)
r.engine2.scale = Math.pow(W_eng/3000, 2/3)
r.engine1.y = 6*lam
r.engine2.y = 6*lam
r.wingrect.scaley = 1-lam
r.wingrect.y = -6 + 5*lam
r.wingtaper.scaley = lam
r.wingtaper.y = 5*lam
"""

def ractivefn(gp):
    sol = gp.solution
    live = "<script>" + ractorfn(sol) + constraints + "</script>"
    display(HTML(live))
    # if you enable the line below, you can try navigating the sliders by sensitivities
    # print sol.table(["cost", "sensitivities"])

with open("flightconditions.gpkit", 'r') as file:
    display(HTML(file.read()))
display(HTML("<style> #ractivecontainer"
             "{position:absolute; height: 0;"
             "right: 0; top: -6em;} </style>"))

gpkit.interactive.widget(gp, ractivefn,
         {"V_{stall,max}": (20, 50, 1),
          "R_{min}": (1e6, 1e7, 0.5e6)})

gpkit.interactive.jswidget(gp, ractorfn,
                           constraints,
           {"V_{stall,max}": (20, 50, 3),
            "R_{min}": (1e6, 1e7, 1e6)})

This concludes the Box example. Try playing around with the sliders up above until you’re bored of this system; then check out one of the other examples. Thanks for reading!

Import CSS for nbviewer¶

If you have a local iPython stylesheet installed, this will add it to the iPython Notebook:

from IPython import utils
from IPython.core.display import HTML
import os
def css_styling():
    """Load default custom.css file from ipython profile"""
    base = utils.path.get_ipython_dir()
    styles = "<style>\n%s\n</style>" % (open(os.path.join(base,'profile_default/static/custom/custom.css'),'r').read())
    return HTML(styles)
css_styling()

Glossary¶

For an alphabetical listing of all commands, check out the Index

The GPkit Package¶

Lightweight GP Modeling Package

For examples please see the examples folder.

Requirements¶

numpy MOSEK or CVXOPT scipy(optional): for complete sparse matrix support sympy(optional): for latex printing in iPython Notebook

Attributes¶

settings : dict

Contains settings loaded from ./env/settings

gpkit.disableUnits()¶

Disables units support in a particular instance of GPkit.

Posynomials created after this is run are incompatible with those created before. If gpkit is imported multiple times, this needs to be run each time.

gpkit.enableUnits()¶

Enables units support in a particular instance of GPkit.

Posynomials created after this is run are incompatible with those created before. If gpkit is imported multiple times, this needs to be run each time.

Subpackages¶

gpkit.interactive¶

gpkit.interactive.plotting module¶

gpkit.interactive.printing module¶

gpkit.interactive.widget module¶

Submodules¶

gpkit.geometric_program¶

Module for creating GP instances.

Example¶

>>> gp = gpkit.GP(cost, constraints, substitutions)

class gpkit.geometric_program.GP(cost, constraints, substitutions={}, solver=None, options={})¶

Bases: gpkit.model.Model

Holds a model and cost function for passing to solvers.

cost : Constraint

Posynomial to minimize when solving

constraints : list of (lists of) Constraints

Constraints to maintain when solving (MonoEQConstraints will be turned into <= and >= constraints)

substitutions : dict {varname: float or int} (optional)

Substitutions to be applied before solving (including sweeps)

solver : str (optional)

Name of solver to use

options : dict (optional)

Options to pass to solver

>>> gp = gpkit.GP(  # minimize
                    0.5*rho*S*C_D*V**2,
                    [   # subject to
                        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,
                        W_w >= W_w_surf + W_w_strc,
                        C_D >= C_D_fuse + C_D_wpar + C_D_ind
                    ], substitutions)
>>> gp.solve()

solve(solver=None, printing=True, skipfailures=False)¶

Solves a GP and returns the solution.

printing : bool (optional)

If True (default), then prints out solver used and time to solve.

solution : dict

A dictionary containing the optimal values for each free variable.

class gpkit.geometric_program.GPSolutionArray¶

Bases: gpkit.small_classes.DictOfLists

DictofLists extended with posynomial substitution.

senssubinto(p)¶

Returns array of each solution’s sensitivity substituted into p

Returns only scalar values.

subinto(p)¶

Returns numpy array of each solution substituted into p.

table(tables=['cost', 'free_variables', 'constants', 'sensitivities'])¶

gpkit.model¶

Module for creating models.

Currently these are only used for GP instances, but they may be further abstractable.

class gpkit.model.Model¶

Bases: object

Abstract class with substituion, loading / saving, and p_idx/A generation

add_constraints(constraints)¶

load(posytuple, printing=True)¶

rm_constraints(constraints)¶

sub(substitutions, val=None, frombase='last', printing=False, replace=False)¶

gpkit.nomials¶

class gpkit.nomials.Constraint(p1, p2)¶

Bases: gpkit.nomials.Posynomial

TODO: Add docstring

class gpkit.nomials.MonoEQConstraint(p1, p2)¶

Bases: gpkit.nomials.Constraint

TODO: Add docstring

class gpkit.nomials.Monomial(exps=None, cs=1, var_locs=None, allow_negative=False, **descr)¶

Bases: gpkit.nomials.Posynomial

TODO: Add docstring

class gpkit.nomials.Posynomial(exps=None, cs=1, var_locs=None, allow_negative=False, **descr)¶

Bases: object

A representation of a posynomial.

exps: tuple of dicts

Exponent dicts for each monomial term

cs: tuple

Coefficient values for each monomial term

var_locs: dict

mapping from variable name to list of indices of monomial terms that variable appears in

Posynomial (if the input has multiple terms) Monomial (if the input has one term)

descr(descr)¶

sub(substitutions, val=None, allow_negative=False)¶

subcmag(substitutions, val=None)¶

to(arg)¶

class gpkit.nomials.VarKey(k=None, **kwargs)¶

Bases: object

Used in the declaration of Posynomials to correspond to each ‘variable name’.

k : object (usually str)

The variable’s name attribute is derived from str(k).

**kwargs :

Any additional attributes, which become the descr attribute (a dict).

VarKey with the given name and descr.

new_unnamed_id = <method-wrapper 'next' of itertools.count object at 0x7fcf836db758>¶

class gpkit.nomials.Variable(*args, **descr)¶

Bases: gpkit.nomials.Monomial

class gpkit.nomials.VectorVariable¶

Bases: gpkit.posyarray.PosyArray

A described vector of singlet Monomials.

length : int

Length of vector.

*args : list

may contain “name” (Strings)

“value” (Iterable) “units” (Strings + Quantity)

and/or “label” (Strings)

**descr : dict

VarKey description

PosyArray of Monomials, each containing a VarKey with name ‘$name_{i}’, where $name is the vector’s name and i is the VarKey’s index.

gpkit.posyarray¶

Module for creating PosyArray instances.

Example¶

>>> x = gpkit.Monomial('x')
>>> px = gpkit.PosyArray([1, x, x**2])

class gpkit.posyarray.PosyArray¶

Bases: numpy.ndarray

A Numpy array with elementwise inequalities and substitutions.

input_array : array-like

>>> px = gpkit.PosyArray([1, x, x**2])

c¶

outer(other)¶

Returns the array and argument’s outer product.

prod()¶

Returns the product of the array.

sub(subs, val=None, allow_negative=False)¶

Substitutes into the array

sum()¶

Returns the sum of the array.

gpkit.small_classes¶

class gpkit.small_classes.CootMatrix¶

Bases: gpkit.small_classes.CootMatrix

A very simple sparse matrix representation.

append(i, j, x)¶

shape = (None, None)¶

tocoo()¶

Converts to a Scipy sparse coo_matrix

tocsc()¶

tocsr()¶

todense()¶

todia()¶

todok()¶

update_shape()¶

gpkit.small_classes.CootMatrixTuple¶

alias of CootMatrix

class gpkit.small_classes.DictOfLists¶

Bases: dict

A hierarchy of dicionaries, with lists at the bottom.

append(sol)¶

Appends a dict (of dicts) of lists to all held lists.

atindex(i)¶

Indexes into each list independently.

toarray(shape=None)¶

Converts all lists into arrays.

class gpkit.small_classes.HashVector¶

Bases: dict

A simple, sparse, string-indexed immutable vector. Inherits from dict.

The HashVector class supports element-wise arithmetic: any undeclared variables are assumed to have a value of zero.

arg : iterable

>>> x = gpkit.nomials.Monomial('x')
>>> exp = gpkit.small_classes.HashVector({x: 2})

class gpkit.small_classes.PosyTuple¶

Bases: tuple

PosyTuple(exps, cs, var_locs, substitutions)

cs¶

Alias for field number 1

exps¶

Alias for field number 0

substitutions¶

Alias for field number 3

var_locs¶

Alias for field number 2

gpkit.small_classes.append_dict(i, o)¶

Recursviely travels dict o and appends items found in i.

gpkit.small_classes.enlist_dict(i, o)¶

Recursviely copies dict i into o, placing non-dict items into lists.

gpkit.small_classes.enray_dict(i, o)¶

Recursively turns lists into numpy arrays.

gpkit.small_classes.index_dict(idx, i, o)¶

Recursviely travels dict i, placing items at idx into dict o.

gpkit.small_scripts¶

gpkit.small_scripts.flatten(ible, classes)¶

Flatten an iterable that contains other iterables

l : Iterable

Top-level container

out : list

List of all objects found in the nested iterables

TypeError

If an object is found whose class was not in classes

gpkit.small_scripts.invalid_types_for_oper(oper, a, b)¶

Raises TypeError for unsupported operations.

gpkit.small_scripts.is_sweepvar(sub)¶

Determines if a given substitution indicates a sweep.

gpkit.small_scripts.latex_num(c)¶

gpkit.small_scripts.locate_vars(exps)¶

From exponents form a dictionary of which monomials each variable is in.

gpkit.small_scripts.mag(c)¶

Return magnitude of a Number or Quantity

gpkit.small_scripts.results_table(data, title, senss=False)¶

gpkit.small_scripts.sort_and_simplify(exps, cs)¶

Reduces the number of monomials, and casts them to a sorted form.

gpkit.small_scripts.unitstr(v, into='%s', options='~')¶

gpkit.substitution¶

Module containing the substitution function

gpkit.substitution.substitution(var_locs, exps, cs, substitutions, val=None)¶

Efficient substituton into a list of monomials.

var_locs : dict

Dictionary of monomial indexes for each variable.

exps : dict

Dictionary of variable exponents for each monomial.

cs : list

Coefficients each monomial.

substitutions : dict

Substitutions to apply to the above.

val : number (optional)

Used to substitute singlet variables.

var_locs_ : dict

Dictionary of monomial indexes for each variable.

exps_ : dict

Dictionary of variable exponents for each monomial.

cs_ : list

Coefficients each monomial.

subs_ : dict

Substitutions to apply to the above.

gpkit.substitution.vectorsub(subs, var, sub, varset)¶

Vectorized substitution via vecmons and Variables.

Citing GPkit¶

If you use GPkit, please cite it using the following bibtex:

@Misc{gpkit,
author = {MIT Convex Optimization Group},
title = {GPkit},
howpublished = {\url{http://github.com/convexopt/gpkit}},
year = {2015},
note = {Version 0.1}
}