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

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:

Unlike most non-linear optimization problems, large GPs can be solved extremely quickly.

If there exists an optimal solution to a GP, it is guaranteed to be globally optimal.

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 at mit.edu.

Mac OS X¶

1. Install Python and build dependencies¶

Install the Python 2.7 version of Anaconda. - Check that Anaconda is installed: in a Terminal window, run python and check that the version string it prints while starting includes “Anaconda”.

If it does not, check that the Anaconda location in .profile in your home directory (you can run vim ~/.profile to read it) corresponds to the location of your Anaconda folder; if it doesn’t, move the Anaconda folder there, and check again in the python startup header.

If you don’t want to install Anaconda, you’ll need gcc, pip, numpy, and scipy, and may find iPython Notebook useful as a modeling environment.

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.

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/

Run pip install ctypesgen --pre in the Terminal (gpkit uses ctypesgen to interface with the MOSEK C bindings)

3. Install GPkit¶

Run pip install gpkit at the command line.

Run pip install ipywidgets for interactive control of models (optional)

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/

Run pip install ctypesgen --pre (gpkit uses ctypesgen to interface with the MOSEK C bindings)

2. Install GPkit¶

_Optional:_ to install gpkit into an isolated python environment, install virtualenv, run virtualenv $DESTINATION_DIR then activate it with source $DESTINATION_DIR/bin/activate.

Run pip install gpkit at the command line.

Run pip install ipywidgets for interactive control of models (optional)

Run python -c "import gpkit.tests; gpkit.tests.run()"

If you want units support, install pint with pip install pint.

You may find iPython Notebook to be useful modeling environment.

Windows¶

1. Install Python dependencies¶

Install the Python 2.7 version of Python (x,y).

Python (x,y) recommends removing any previous installations of Python before installation.

Make sure to check the cvxopt checkbox under “Choose components” during installation.

2. (optional) Install the MOSEK GP solver¶

CVXOPT is included with Python (x,y) and does not need to be installed

Installing CVXOPT with Anaconda or another Python distribution can be difficult, which is why we reccomend Python (x,y).

Download MOSEK, then:

Follow these steps for Windows.

Request an academic license file and put it in ~/mosek/

To use the MOSEK C bindings solver: - Make sure “gcc” is on your system path (that is, you can type gcc into a command prompt and not get “executable not found”) - Run pip install ctypesgen --pre in the Command Prompt (gpkit uses ctypesgen to interface with the MOSEK C bindings)

3. Install GPkit¶

Run pip install gpkit at an Anaconda Command Prompt.

Run pip install ipywidgets for interactive control of models (optional)

Run python -c "import gpkit.tests; gpkit.tests.run()"

If attempting to run the tests results in ValueError: Unknown solver ''. and python -c "import gpkit" prints “Could not find settings file”, then run python -c "import gpkit.build; gpkit.build.build_gpkit()", to look for and install solvers. After doing so (it should say which solvers have been found), run python -c "import gpkit.tests; gpkit.tests.run()" again.

If you want units support, install pint with pip install pint.

Updating GPkit between releases¶

Active developers may wish to install the latest GPkit <http://github.com/hoburg/gpkit> directly from the source code on Github. To do so,

Run pip uninstall gpkit to uninstall your existing GPkit.

Run git clone https://github.com/hoburg/gpkit.git to clone the GPkit repository, or cd gpkit; git pull origin master; cd .. to update your existing repository.

Run pip install -e gpkit to reinstall GPkit.

Run python -c "import gpkit.tests; gpkit.tests.run()" to test your installation.
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:

Open your command line interface (terminal/Command Prompt).

Open a Python interpreter. This can be done by typing python (or ipython if installed).

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, Model
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.
# NOTE: starting a Python string with 'r' makes the backslashes literal,
#       which is useful for LaTeX strings.
rho = Variable(r'\rho', 1.225, 'kg/m^3', 'Density of air at sea level')
In the example above, the key name r'\rho' is for LaTeX printing (described later). The unit and description arguments are optional.
#Declare pi equal to 3.14
pi = Variable(r'\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 Model-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 Model¶

The Model class represents an optimization 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 = Model(objective, constraints)
Solving the Model¶
sol = gp.solve()
Printing Results¶
print sol.table()
print "The x dimension is %s." % (sol(x))
Advanced Commands¶
Feasibility Analysis¶

If your Model doesn’t solve, you can automatically find the nearest feasible version of it with the Model.feasibility() command, as shown below. The feasible version can either involve relaxing all constraints by the smallest number possible (that is, dividing the less-than side of every constraint by the same number), relaxing each constraint by its own number and minimizing the product of those numbers, or changing each constant by the smallest total percentage possible.
from gpkit import Variable, Model, PosyArray
x = Variable("x")
x_min = Variable("x_min", 2)
x_max = Variable("x_max", 1)
m = Model(x, [x <= x_max, x >= x_min])
# m.solve()  # raises a RuntimeWarning!
feas = m.feasibility()

# USING OVERALL
m.constraints = PosyArray(m.signomials)/feas["overall"]
m.solve()

# USING CONSTRAINTS
m = Model(x, [x <= x_max, x >= x_min])
m.constraints = PosyArray(m.signomials)/feas["constraints"]
m.solve()

# USING CONSTANTS
m = Model(x, [x <= x_max, x >= x_min])
m.substitutions.update(feas["constants"])
m.solve()
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 SolutionArray’s sens() method, as in this example:
import gpkit
x = gpkit.Variable("x")
x_min = gpkit.Variable("x_{min}", 2)
sol = gpkit.Model(x, [x_min <= x]).solve()
assert sol.sens(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.Model(x, [x_squared_min <= x**2]).solve()
assert sol.sens(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(m)
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 Model 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 Model’s cost and constraints before the substitutions argument is (and hence values are supplanted by any later substitutions).

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(solution["variables"]).

You can also substitute by just calling the solution, i.e. 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 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 Model 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 parallelized.

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, Model

x = Variable("x")
y = VectorVariable(2, "y")
m = Model(x, [x >= y.prod()])
m.substitutions.update({y: ('sweep', [[2, 3], [5, 7, 11]])})
a = m.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.tools.composite_objective(L+W, W**-1 * L**-3,
                                      normsub={L:10, W: 10},
                                      sweep=co_sweep)

m = gpkit.Model(obj, eqns)
m.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:

maximize the number of variables in $g$,

prioritize variables that are in the objective,

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
with gpkit.SignomialsEnabled():
    constraints = [x >= 1-y, y <= 0.1]

# create and solve the SP
m = gpkit.Model(x, constraints)
sol = m.localsolve(verbosity=1)
assert abs(sol(x) - 0.9) < 1e-6
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, 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
from gpkit.shortcuts import *


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)
L, EI, q = sol("L"), sol("EI"), sol("q")
N = len(q)
q = q[0]  # assume uniform loading for the check below
x = np.linspace(0, L, N)  # position along beam
w_gp = sol("w")  # deflection along beam
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

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.
Glossary¶

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

The GPkit Package¶

Subpackages¶

gpkit.interactive¶

gpkit.interactive.plotting module¶

gpkit.interactive.pretty_print module¶

gpkit.interactive.printing module¶

gpkit.interactive.ractor module¶

gpkit.interactive.widget module¶

Submodules¶

gpkit.model¶

gpkit.solution_array module¶

gpkit.geometric_program¶

gpkit.signomial_program¶

gpkit.nomials¶

gpkit.posyarray¶

gpkit.small_classes¶

gpkit.small_scripts¶

gpkit.substitution¶

gpkit.tools module¶

gpkit.variables module¶

gpkit.varkey module¶
Citing GPkit¶

If you use GPkit, please cite it with the following bibtex:
@Misc{gpkit,
      author={Edward Burnell and Warren Hoburg},
      title={GPkit},
      howpublished={\url{https://github.com/hoburg/gpkit}},
      year={2015},
      note={Version 0.3.4}
     }
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.3¶

Integrated GP and SP creation under the Model class

Improved and simplified under-the-hood internals of GPs and SPs

New experimental SP heuristic

Improved test coverage

Handles vectors which are partially constants, partially free

Simplified interaction with Model objects and made it more pythonic

Added SP “step” method to allow single-stepping through an SP

Isolated and corrected some solver-specific behavior

Fully allowed substitutions of variables for 0 (commit 4631255)

Use “with” to create a signomials environment (commit cd8d581)

Continuous integration improvements, thanks @galbramc !

Not counting subpackages, went from 2200 to 2400 lines of code (additions were mostly longer error messages) and from 650 to 1050 lines of docstrings and comments.

Add automatic feasibility-analysis methods to Model and GP

Simplified solver logging and printing, making it easier to access solver output.

Version 0.2¶

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)