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 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 decision variables, and \(a_{1..n}\) are real exponents. For example, taking \(x\), \(y\) and \(z\) to be positive variables, the expressions

\[7x \qquad 4xy^2z \qquad \frac{2x}{y^2z^{0.3}} \qquad \sqrt{2xy}\]

are all monomials. 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}\]

For example, the expressions

\[x^2 + 2xy + 1 \qquad 7xy + 0.4(yz)^{-1/3} \qquad 0.56 + \frac{x^{0.7}}{yz}\]

are all posynomials. Alternatively, monomials can be defined as the subset of posynomials having only one term. Using \(f_i\) to represent a monomial and \(g_i\) to represent a posynomial, a GP in standard form is written as:

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

Boyd et. al. give the following example of a GP in standard form:

\[\begin{split}\begin{array}[llll]\text{} \text{minimize} & x^{-1}y^{-1/2}z^{-1} + 2.3xz + 4xyz \\ \text{subject to} & (1/3)x^{-2}y^{-2} + (4/3)y^{1/2}z^{-1} \leq 1 \\ & x + 2y + 3z \leq 1 \\ & (1/2)xy = 1 \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.

What are Signomials / Signomial Programs?¶

When the coefficients in a posynomial are allowed to be negative (but the variables stay strictly positive), that is called a Signomial. A Signomial Program has signomial constraints, and while they cannot be solved as quickly or to global optima, because they build on the structure of a GP they can be solved more quickly than a generic nonlinear program. More information on Signomial Programs can be found under “Advanced Commands”.

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.

GPkit Overview¶

GPkit is a Python package for defining and manipulating geometric programming (GP) models, abstracting away the backend solver.

The primary goal of GPkit is to make it easy to create share and explore Geometric Programming models. Being fast and mathematically correct tend to align well with this goal.

Symbolic expressions¶

GPkit is a limited symbolic algebra language, allowing only for the creation of Geometric Program compatible equations (or Signomial Program compatible ones, if signomial programming is enabled). As mentioned in “Geometric Programming 101”, one can look at monomials as posynomials with a single term, and posynomials as signomials that have only positive coefficients. The inheritance structure of these objects in GPkit follows this mathematical basis.

Substitution¶

The Varkey object in the graph above is not a algebraic expression, but what GPkit uses as a variable’s “name”. It carries the LaTeX representation of a variable and its units, as well as any other information the user wishes to associate with a variable. The use of VarKeys as opposed to numeric indexing is an important part of the GPkit framework, because it allows a user to keep their model entirely symbolic until they wish to solve it.

This means that any expression or Model object can replace any instance of a variable (as represented by a VarKey) with a number, new VarKey, or even an entire Monomial at any time with the .sub() method.

Model objects¶

In GPkit, a Model object represents a symbolic problem declaration. That problem may be either GP-compatible or SP-compatible. To avoid confusion, calling the solve() method on a model will either attempt to solve it for a global optimum (if it’s a GP) or return an error immediately (if it’s an SP). Similarly, calling localsolve() will either start the process of SP-solving (stepping through a sequence of GP-approximations) or return an error for GP-compatible Models. This framework is illustrated below.

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

Advanced Commands¶

Sensitivities and dual variables¶

When a GP is solved, the solver returns not just the optimal value for the problem’s variables (known as the “primal solution”) but also, as a side effect of the solving process, the effect that scaling the less-than side of each constraint would have on the overall objective (called the “dual solution”, “shadow prices”, or “posynomial sensitivities”).

Using variable sensitivities¶

GPkit uses this dual solution to compute the sensitivities of each variable, which can be accessed most easily using a GPSolutionArray’s senssubinto() method, as in this example:

import gpkit
x = gpkit.Variable("x")
x_min = gpkit.Variable("x_{min}", 2)
sol = gpkit.GP(x, [x_min <= x]).solve()
assert sol.senssubinto(x_min) == 1

These sensitivities are actually log derivatives (\(\frac{d \mathrm{log}(y)}{d \mathrm{log}(x)}\)); whereas a regular derivative is a tangent line, these are tangent monomials, so the 1 above indicates that x_min has a linear relation with the objective. This is confirmed by a further example:

import gpkit
x = gpkit.Variable("x")
x_squared_min = gpkit.Variable("x^2_{min}", 2)
sol = gpkit.GP(x, [x_squared_min <= x**2]).solve()
assert sol.senssubinto(x_squared_min) == 2

Plotting variable sensitivities¶

Sensitivities are a useful way to evaluate the tradeoffs in your model, as well as what aspects of the model are driving the solution and should be examined. To help with this, GPkit has an automatic sensitivity plotting function that can be accessed as follows:

from gpkit.interactive.plotting import sensitivity_plot
sensitivity_plot(gp)

Which produces the following plot:

In this plot, steep lines that go up to the right are variables whose increase sharply increases (makes worse) the objective. Steep lines going down to the right are variables whose increase sharply decreases (improves) the objective.

Substitutions¶

Substitutions are a general-purpose way to change every instance of one variable into either a number or another variable.

Substituting into Posynomials, PosyArrays, and GPs¶

The examples below all use Posynomials and PosyArrays, but the syntax is identical for GPs (except when it comes to sweep variables).

# adapted from t_sub.py / t_NomialSubs / test_Basic
from gpkit import Variable
x = Variable("x")
p = x**2
assert p.sub(x, 3) == 9
assert p.sub(x.varkeys["x"], 3) == 9
assert p.sub("x", 3) == 9

Here the variable x is being replaced with 3 in three ways: first by substituting for x directly, then by substituting for the VarKey("x"), then by substituting the string “x”. In all cases the substitution is understood as being with the VarKey: when a variable is passed in the VarKey is pulled out of it, and when a string is passed in it is used as an argument to the Posynomial’s varkeys dictionary.

Substituting multiple values¶

# adapted from t_sub.py / t_NomialSubs / test_Vector
from gpkit import Variable, VectorVariable
x = Variable("x")
y = Variable("y")
z = VectorVariable(2, "z")
p = x*y*z
assert all(p.sub({x: 1, "y": 2}) == 2*z)
assert all(p.sub({x: 1, y: 2, "z": [1, 2]}) == z.sub(z, [2, 4]))

To substitute in multiple variables, pass them in as a dictionary where the keys are what will be replaced and values are what it will be replaced with. Note that you can also substitute for VectorVariables by their name or by their PosyArray.

Substituting with nonnumeric values¶

You can also substitute in sweep variables (see Sweeps), strings, and monomials:

# adapted from t_sub.py / t_NomialSubs
from gpkit import Variable
from gpkit.small_scripts import mag

x = Variable("x", "m")
xvk = x.varkeys.values()[0]
descr_before = x.exp.keys()[0].descr
y = Variable("y", "km")
yvk = y.varkeys.values()[0]
for x_ in ["x", xvk, x]:
    for y_ in ["y", yvk, y]:
        if not isinstance(y_, str) and type(xvk.units) != str:
            expected = 0.001
        else:
            expected = 1.0
        assert abs(expected - mag(x.sub(x_, y_).c)) < 1e-6
if type(xvk.units) != str:
    # this means units are enabled
    z = Variable("z", "s")
    # y.sub(y, z) will raise ValueError due to unit mismatch

Note that units are preserved, and that the value can be either a string (in which case it just renames the variable), a varkey (in which case it changes its description, including the name) or a Monomial (in which case it substitutes for the variable with a new monomial).

Substituting with replacement¶

Any of the substitutions above can be run with p.sub(*args, replace=True) to clobber any previously-substitued values.

Fixed Variables¶

When a GP is created, any fixed Variables are used to form a dictionary: {var: var.descr["value"] for var in self.varlocs if "value" in var.descr}. This dictionary in then substituted into the GP’s cost and constraints before the substitutions argument.

Substituting from a GP solution array¶

gp.solution.subinto(p) will substitute the solution(s) for variables into the posynomial p, returning a PosyArray. For a non-swept solution, this is equivalent to p.sub(gp.solution["variables"]).

You can also substitute by just calling the solution, i.e. gp.solution(p). This returns a numpy array of just the coefficients (c) of the posynomial after substitution, and will raise a` ValueError` if some of the variables in p were not found in gp.solution.

Sweeps¶

Declaring Sweeps¶

Sweeps are useful for analyzing tradeoff surfaces. A sweep “value” is an Iterable of numbers, e.g. [1, 2, 3]. Variables are swept when their substitution value takes the form ('sweep', Iterable), (e.g. 'sweep', np.linspace(1e6, 1e7, 100)). During variable declaration, giving an Iterable value for a Variable is assumed to be giving it a sweeep value: for example, x = Variable("x", [1, 2, 3]. Sweeps can also be declared during later substitution (gp.sub("x", ('sweep', [1, 2, 3])), or if the variable was already substituted for a constant, gp.sub("x", ('sweep', [1, 2, 3]), replace=True)).

Solving Sweeps¶

A GP with sweeps will solve for all possible combinations: e.g., if there’s a variable x with value ('sweep', [1, 3]) and a variable y with value ('sweep', [14, 17]) then the gp will be solved four times, for \((x,y)\in\left\{(1, 14),\ (1, 17),\ (3, 14),\ (3, 17)\right\}\). The returned solutions will be a one-dimensional array (or 2-D for vector variables), accessed in the usual way. Sweeping Vector Variables

Vector variables may also be substituted for: y = VectorVariable(3, "y", value=('sweep' ,[[1, 2], [1, 2], [1, 2]]) will sweep \(y\ \forall~y_i\in\left\{1,2\right\}\).

Parallel Sweeps¶

During a normal sweep, each result is independent, so they can be run in parallel. To use this feature, run $ ipcluster start at a terminal: it will automatically start a number of iPython parallel computing engines equal to the number of cores on your machine, and when you next import gpkit you should see a note like Using parallel execution of sweeps on 4 clients. If you do, then all sweeps performed with that import of gpkit will be parellelized.

This parallelization sets the stage for gpkit solves to be outsourced to a server, which may be valuable for faster results; alternately, it could allow the use of gpkit without installing a solver.

Linked Sweeps¶

Some constants may be “linked” to another sweep variable. This can be represented by a Variable whose value is ('sweep', fn), where the arguments of the function fn are stored in the Varkeys’s args attribute. If you declare a variables value to be a function, then it will assume you meant that as a sweep value: for example, a_ = gpkit.Variable("a_", lambda a: 1-a, "-", args=[a]) will create a constant whose value is always 1 minus the value of a (valid for values of a less than 1). Note that this declaration requires the variable a to already have been declared.

Example Usage¶

# code from t_GPSubs.test_VectorSweep in tests/t_sub.py
from gpkit import Variable, VectorVariable, GP

x = Variable("x")
y = VectorVariable(2, "y")
gp = GP(x, [x >= y.prod()])
gp.sub(y, ('sweep', [[2, 3], [5, 7, 11]]))
a = gp.solve(printing=False)["cost"]
b = [10, 14, 22, 15, 21, 33]
assert all(abs(a-b)/(a+b) < 1e-7)

Composite Objectives¶

Given \(n\) posynomial objectives \(g_i\), you can sweep out the problem’s Pareto frontier with the composite objective:

\(g_0 w_0 \prod_{i\not=0} v_i + g_1 w_1 \prod_{i\not=1} v_i + ... + g_n \prod_i v_i\)

where \(i \in 0 ... n-1\) and \(v_i = 1- w_i\) and \(w_i \in [0, 1]\)

GPkit has the helper function composite_objective for constructing these.

Example Usage¶

import numpy as np
import gpkit

L, W = gpkit.Variable("L"), gpkit.Variable("W")

eqns = [L >= 1, W >= 1, L*W == 10]

co_sweep = [0] + np.logspace(-6, 0, 10).tolist()

obj = gpkit.composite_objective(L+W, W**-1 * L**-3,
                                normsub={L:10, W: 10},
                                sweep=co_sweep)

gp = gpkit.GP(obj, eqns)
gp.solve()

The normsub argument specifies an expected value for your solution to normalize the different \(g_i\) (you can also do this by hand). The feasibility of the problem should not depend on the normalization, but the spacing of the sweep will.

The sweep argument specifies what points between 0 and 1 you wish to sample the weights at. If you want different resolutions or spacings for different weights, the sweeps argument accepts a list of sweep arrays.

Signomial Programming¶

Signomial programming finds a local solution to a problem of the form:

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

where each \(f\) is monomial while each \(g\) and \(h\) is a posynomial.

This requires multiple solutions of geometric programs, and so will take longer to solve than an equivalent geometric programming formulation.

The specification of the signomial problem affects its solve time in a nuanced way: gpkit.SP(x, [x >= 0.1, x+y >= 1, y <= 0.1]).localsolve() takes about four times as many iterations to solve as gpkit.SP(x, [x >= 1-y, y <= 0.1]).localsolve(), despite the two formulations being arithmetically equivalent.

In general, when given the choice of which variables to include in the positive-posynomial / \(g\) side of the constraint, the modeler should:

1. maximize the number of variables in \(g\),
2. prioritize variables that are in the objective,
3. then prioritize variables that are present in other constraints.

The syntax SP.localsolve is chosen to emphasize that signomial programming returns a local optimum. For the same reason, calling SP.solve will raise an error.

By default, signomial programs are first solved conservatively (by assuming each \(h\) is equal only to its constant portion) and then become less conservative on each iteration.

Example Usage¶

"""Adapted from t_SP in tests/t_geometric_program.py"""
import gpkit

# Decision variables
x = gpkit.Variable('x')
y = gpkit.Variable('y')

# must enable signomials for subtraction
gpkit.enable_signomials()

# create and solve the SP
sp = gpkit.SP(x, [x >= 1-y, y <= 0.1])
sol = sp.localsolve(printing=False)
assert abs(sol(x) - 0.9) < 1e-6

gpkit.disable_signomials()

When using the localsolve method, the reltol argument specifies the relative tolerance of the solver: that is, by what percent does the solution have to improve between iterations? If any iteration improves less than that amount, the solver stops and returns its value.

If you wish to start the local optimization at a particular point \(x_k\), however, you may do so by putting that position (a dictionary formatted as you would a substitution) as the xk argument.

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, 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(printing=False)

# 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, 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(printing=False)

# Print results table
print sol.table()

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 

Swept variables
---------------

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

Constant and swept variable 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, 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.table()

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        

Swept variables
---------------

Constants
---------
   M : 100   [kg]      Mass of Water in the Tank   
\rho : 1000  [kg/m**3] Density of Water in the Tank

Constant and swept variable 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.

from gpkit.shortcuts import *
import numpy as np

# Define the constants in the problem
CDA0 = Var('(CDA_0)', 0.0306, 'm^2', label='Fuselage drag area')
C_Lmax = Var('C_{L,max}', 2.0, label='Maximum C_L, flaps down')
e = Var('e', 0.96, label='Oswald efficiency factor')
k = Var('k', 1.2, label='Form factor')
mu = Var('\\mu', 1.78E-5, 'kg/m/s', label='Viscosity of air')
N_lift = Var('N_{lift}', 2.5, label='Ultimate load factor')
rho = Var('\\rho', 1.23, 'kg/m^3', label='Density of air')
Sw_S = Var('\\frac{S_{wet}}{S}', 2.05, label='Wetted area ratio')
tau = Var('\\tau', 0.12, label='Airfoil thickness-to-chord ratio')
V_min = Var('V_{min}', 22, 'm/s', label='Desired landing speed')
W_0 = Var('W_0', 4940, 'N', label='Aircraft weight excluding wing')
cww1 = Var('cww1', 45.42, 'N/m^2', 'Wing weight area factor')
cww2 = Var('cww2', 8.71E-5, '1/m', 'Wing weight bending factor')

# Define decision variables
A = Var('A', label='Aspect ratio')
C_D = Var('C_D', label='Drag coefficient')
C_L = Var('C_L', label='Lift coefficient')
C_f = Var('C_f', label='Skin friction coefficient')
S = Var('S', 'm^2', label='Wing planform area')
Re = Var('Re', label='Reynolds number')
W = Var('Re', 'N', label='Total aircraft weight')
W_w = Var('W_w', 'N', label='Wing weight')
V = Var('V', 'm/s', label='Cruise velocity')

# Define objective function
objective = 0.5 * rho * V**2 * C_D * S

# Define constraints
constraints = [C_f * Re**0.2 >= 0.074,
               C_D >= CDA0 / S + k * C_f * Sw_S + C_L**2 / (np.pi * A * e),
               0.5 * rho * V**2 * C_L * S >= W,
               W >= W_0 + W_w,
               W_w >= (cww1*S +
                       cww2 * N_lift * A**1.5 * (W_0 * W * S)**0.5 / tau),
               0.5 * rho * V_min**2 * C_Lmax * S >= W,
               Re == rho * V * (S/A)**0.5 / mu]

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

# Solve the GP
sol = gp.solve()

# Print the solution table
print sol.table()

The output is

Using solver 'cvxopt'
Solving for 9 variables.
Solving took 0.0629 seconds.

Cost
----
 255 [kg*m/s**2] 

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 : 7189       [N]    Total aircraft weight    
 Re : 2.598e+06         Reynolds number          
  S : 12.08      [m**2] Wing planform area       
  V : 38.55      [m/s]  Cruise velocity          
W_w : 2249       [N]    Wing weight              

Swept variables
---------------

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                     

Constant and swept variable 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                     

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.disable_signomials()¶

Disables signomial support in a particular instance of GPkit.

gpkit.disable_units()¶

Disables units support in a particular instance of GPkit.

Posynomials created after calling this are incompatible with those created before.

If gpkit is imported multiple times, this needs to be run each time.

The correct way to call this is:

import gpkit gpkit.disable_units()

The following will not have the intended effect:

from gpkit import disable_units disable_units()

gpkit.enable_signomials()¶

Enables signomial support in a particular instance of GPkit.

gpkit.enable_units()¶

Enables units support in a particular instance of GPkit.

Posynomials created after calling this 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¶

gpkit.interactive.ractor module¶

Submodules¶

gpkit.geometric_program¶

Module for creating GeometricProgram instances.

Example¶

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

class gpkit.geometric_program.GPSolutionArray¶

Bases: gpkit.small_classes.DictOfLists

DictofLists extended with posynomial substitution.

getvars(*args)¶

senssubinto(p)¶

Returns array of each solution’s sensitivity substituted into p

Returns only scalar values.

subinto(p)¶

Returns PosyArray of each solution substituted into p.

table(tables=['cost', 'freevariables', 'sweptvariables', 'constants', 'sensitivities'], fixedcols=True)¶

class gpkit.geometric_program.GeometricProgram(*args, **kwargs)¶

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.GeometricProgram(  # 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()

model_nums = defaultdict(<type 'int'>, {})¶

solve(*args, **kwargs)¶

test(solver=None, printing=False, skipfailures=False)¶

vars(*args)¶

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

checkbounds()¶

genA(allpos=True, printing=True)¶

load(posytuple)¶

rm_constraints(constraints)¶

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

Substitutes into a model, modifying it in place.

gp.sub({‘x’: 1, y: 2}) gp.sub(x, 3, replace=True)

substitutions : dict or key

Either a dictionary whose keys are strings, Variables, or VarKeys, and whose values are numbers, or a string, Variable or Varkey.

val : number (optional)

If the substitutions entry is a single key, val holds the value

frombase : string (optional)

Which model state to update. By default updates the current state.

replace : boolean (optional, default is True)

Whether the substitution should only subtitute currently unsubstituted values (False) or should also make replacements of current substitutions (True).

No return value: the model is modified in place.

gpkit.nomials¶

Signomial, Posynomial, Monomial, Constraint, & MonoEQCOnstraint classes

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

Bases: gpkit.nomials.Posynomial

A constraint of the general form monomial > posynomial Stored internally (exps, cs) as a single Posynomial (1 >= self) Usually initialized via operator overloading, e.g. cc = y**2 >= 1 + x Additionally stores input format (lhs vs rhs) in self.right and self.left Form is self.left >= self.right.

TODO: this documentation needs to address Signomial Constraints.

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

Bases: gpkit.nomials.Constraint

A Constraint of the form Monomial == Monomial. Stored internally as a Monomial Constraint, (1 == self).

class gpkit.nomials.Monomial(exps=None, cs=1, varlocsandkeys=None, require_positive=True, **descr)¶

Bases: gpkit.nomials.Posynomial

A Posynomial with only one term

Same as Signomial. Note: Monomial historically supported several different init formats

These will be deprecated in the future, replaced with a single __init__ syntax, same as Signomial.

class gpkit.nomials.Posynomial(exps=None, cs=1, varlocsandkeys=None, require_positive=True, **descr)¶

Bases: gpkit.nomials.Signomial

A Signomial with strictly positive cs

Same as Signomial. Note: Posynomial historically supported several different init formats

These will be deprecated in the future, replaced with a single __init__ syntax, same as Signomial.

class gpkit.nomials.Signomial(exps=None, cs=1, varlocsandkeys=None, require_positive=True, **descr)¶

Bases: object

A representation of a signomial.

exps: tuple of dicts

Exponent dicts for each monomial term

cs: tuple

Coefficient values for each monomial term

varlocsandkeys: dict

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

require_positive: bool

If True and signomials not enabled, c <= 0 will raise ValueError

Signomial Posynomial (if the input has only positive cs) Monomial (if the input has one term and only positive cs)

descr(descr)¶

diff(wrt)¶

Derivative of this with respect to a Variable

wrt (Variable):

Variable to take derivative with respect to

Signomial (or Posynomial or Monomial)

mono_approximation(x0)¶

prod()¶

sub(substitutions, val=None, require_positive=True)¶

Returns a nomial with substitued values.

3 == (x**2 + y).sub({‘x’: 1, y: 2}) 3 == (x).gp.sub(x, 3)

substitutions : dict or key

Either a dictionary whose keys are strings, Variables, or VarKeys, and whose values are numbers, or a string, Variable or Varkey.

val : number (optional)

If the substitutions entry is a single key, val holds the value

require_positive : boolean (optional, default is True)

Controls whether the returned value can be a signomial.

Returns substituted nomial.

subcmag(substitutions, val=None)¶

sum()¶

to(arg)¶

value¶

Self, with values substituted for variables that have values

float, if no symbolic variables remain after substitution (Monomial, Posynomial, or Signomial), otherwise.

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¶

left¶

Self, sampled one index down, with zeropad

outer(other)¶

Returns the array and argument’s outer product.

right¶

Self, sampled one index up, with zeropad

sub(subs, val=None, require_positive=True)¶

Substitutes into the array

gpkit.small_classes¶

Miscellaneous small classes

class gpkit.small_classes.CootMatrix¶

Bases: gpkit.small_classes.CootMatrix

A very simple sparse matrix representation.

append(row, col, data)¶

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(*args, **kwargs)¶

Bases: dict

A simple, sparse, string-indexed 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(exps, cs, varlocs, substitutions)¶

Bases: tuple

cs¶

Alias for field number 1

exps¶

Alias for field number 0

substitutions¶

Alias for field number 3

varlocs¶

Alias for field number 2

gpkit.small_classes.append_dict(i, o)¶

Recursviely travels dict o and appends items found in i.

class gpkit.small_classes.count¶

Bases: object

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.diff(p, vk)¶

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.isequal(a, b)¶

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.mono_approx(p, x0)¶

gpkit.small_scripts.results_table(data, title, minval=0, printunits=True, fixedcols=True, varfmt='%s : ', valfmt='%-.4g ', vecfmt='%-8.3g')¶

Pretty string representation of a dict of VarKeys Iterable values are handled specially (partial printing)

data: dict whose keys are VarKey’s

data to represent in table

title: string minval: float

skip values with all(abs(value)) < minval

printunits: bool fixedcols: bool

if True, print rhs (val, units, label) in fixed-width cols

varfmt: string

format for variable names

valfmt: string

format for scalar values

vecfmt: string

format for vector values

gpkit.small_scripts.sort_and_simplify(exps, cs)¶

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

gpkit.small_scripts.unitstr(units, into='%s', options='~', dimless='-')¶

gpkit.substitution¶

Module containing the substitution function

gpkit.substitution.getsubs(varkeys, varlocs, substitutions)¶

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

Efficient substituton into a list of monomials.

varlocs : dict

Dictionary mapping variables to lists of monomial indices.

exps : Iterable of dicts

Dictionary mapping variables to exponents, for each monomial.

cs : list

Coefficient for each monomial.

substitutions : dict

Substitutions to apply to the above.

val : number (optional)

Used to substitute singlet variables.

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

Citing GPkit¶

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

@Misc{gpkit,
      author={MIT Department of Aeronautics and Astronautics},
      title={GPkit},
      howpublished={\url{https://github.com/convexopt/gpkit}},
      year={2015},
      note={Version 0.2}
     }

Acknowledgements¶

We thank the following contributors for helping to improve GPkit:

• Marshall Galbraith for setting up continuous integration.
• Stephen Boyd for inspiration and suggestions.

Release Notes¶

This page lists the changes made in each point version of gpkit.

Version 0.2.0¶

• Various bug fixes
• Python 3 compatibility
• Added signomial programming support (alpha quality, may be wrong)
• Added composite objectives
• Parallelized sweeping
• Better table printing
• Linked sweep variables
• Better error messages
• Closest feasible point capability
• Improved install process (no longer requires ctypesgen; auto-detects MOSEK version)
• Added examples: wind turbine, modular GP, examples from 1967 book, maintenance (part replacement)
• Documentation grew by ~70%
• Added Advanced Commands section to documentation
• Many additional unit tests (more than doubled testing lines of code)