GPkit Documentation¶
Welcome!
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 a GP, a monomial is defined as:
\[f(x) = c x_1^{a_1} x_2^{a_2} ... x_n^{a_n}\]where \(c\) is a positive constant, \(x_{1..n}\) are the decision variables, and \(a_{1..n}\) are real exponents.
Building on this, a posynomial is defined as a sum of monomials:
\[g(x) = \sum_{k=1}^K c_k x_1^{a_1k} x_2^{a_2k} ... x_n^{a_nk}\]Using these definitions, a GP in Standard Form is written as:
\[\begin{split}\begin{array}[lll]\text{} \text{minimize} & f_0(x) & \\ \text{subject to} & f_i(x) = 1, & i = 1,....,m \\ & g_i(x) \leq 1, & i = 1,....,n \end{array}\end{split}\]Why are GPs special?¶
Geometric programs have several powerful properties:
- 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.
Where can I learn more?¶
To learn more about GPs, refer to the following resources:
- A tutorial on geometric programming, by S. Boyd, S.J. Kim, L. Vandenberghe, and A. Hassibi.
- Convex optimization, by S. Boyd and L. Vandenberghe.
Installation Instructions¶
If you encounter any bugs during installation, please email gpkit@mit.edu.
Mac OS X¶
1. Install Python and build dependencies¶
- Install the Python 2.7 version of Anaconda.
- If you don’t want to install Anaconda, you’ll need gcc and pip, and may find sympy, scipy, and iPython Notebook useful.
- If which gcc does not return anything, install the Apple Command Line Tools.
- _Optional:_ to install gpkit into an isolated python environment you can create a new conda virtual environment with conda create -n gpkit anaconda and activate it with source activate gpkit.
- Run pip install ctypesgen in the Terminal.
2. Install either the MOSEK or CVXOPT GP solvers¶
- Download CVXOPT, then:
- Read the official instructions and requirements
- In the Terminal, navigate to the cvxopt folder
- Run python setup.py install
- Download MOSEK, then:
- Move the mosek folder to your home directory
- Follow these steps for Mac.
- Request an academic license file and put it in ~/mosek/
3. Install GPkit¶
- Run pip install gpkit at the command line.
- Run python -c "import gpkit.tests; gpkit.tests.run()"
- If you want units support, install pint with pip install pint.
Linux¶
1. Install either the MOSEK or CVXOPT GP solvers¶
- Download CVXOPT, then:
- Read the official instructions and requirements
- In a terminal, navigate to the cvxopt folder
- Run python setup.py install
- Download MOSEK, then:
- Move the mosek folder to your home directory
- Follow these steps for Linux.
- Request an academic license file and put it in ~/mosek/
2. Install GPkit¶
- _Optional:_ to install gpkit into an isolated python environment, install virtualenv, run virtualenv $DESTINATION_DIR then activate it with source activate $DESTINATION_DIR/bin.
- Run pip install ctypesgen at the command line.
- Run pip install gpkit at the command line.
- Run python -c "import gpkit.tests; gpkit.tests.run()"
- If you want units support, install pint with pip install pint.
- You may find sympy, scipy, and iPython Notebook to be useful additional packages as well.
Windows¶
1. Install Python dependencies¶
- Install the Python 2.7 version of Anaconda.
- If you don’t want to install Anaconda, you’ll need gcc and pip, and may find sympy, scipy, and iPython Notebook useful.
- _Optional:_ to install gpkit into an isolated python environment you can create a new conda virtual environment with conda create -n gpkit anaconda and activate it with source activate gpkit.
- Run pip install ctypesgen at an Anaconda Command Prompt.
2. Install either the MOSEK or CVXOPT GP solvers¶
Download CVXOPT, then follow these steps to install a linear algebra library
- Download MOSEK, then:
- Follow these steps for Windows.
- Request an academic license file and put it in ~/mosek/
3. Install GPkit¶
- Run pip install gpkit at an Anaconda Command Prompt.
- Run python -c "import gpkit.tests; gpkit.tests.run()"
- If you want units support, install pint with pip install pint.
Getting Started¶
GPkit is a Python package. We assume basic familiarity with Python. If you are new to Python take a look at Learn Python.
GPkit is also a command line tool. This means that you need to be in the terminal (OS X/Linux) or command prompt (Windows) to use it. If you are not familiar with working in the command line, check out this Learn Code the Hard Way tutorial.
The first thing to do is install GPkit . Once you have done this, you can start using GPkit in 3 easy steps:
- 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.pyAgain, 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, GPDeclaring 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.PosynomialDeclaring 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.ConstraintDeclaring 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()Examples¶
A Trivial GP¶
from gpkit import Variable, GP # Decision variable x = Variable('x') # Constraint constraints = [x >= 1] # Objective (to minimize) objective = x # Formulate the GP gp = GP(objective, constraints) # Solve the GP sol = gp.solve() # Print results table print sol.table()Maximizing the Volume of a Box¶
from gpkit import Variable, GP # Parameters alpha = Variable("alpha", 2, "-", "lower limit, wall aspect ratio") beta = Variable("beta", 10, "-", "upper limit, wall aspect ratio") gamma = Variable("gamma", 2, "-", "lower limit, floor aspect ratio") delta = Variable("delta", 10, "-", "upper limit, floor aspect ratio") A_wall = Variable("A_{wall}", 200, "m^2", "upper limit, wall area") A_floor = Variable("A_{floor}", 50, "m^2", "upper limit, floor area") # Decision variables h = Variable("h", "m", "height") w = Variable("w", "m", "width") d = Variable("d", "m", "depth") #Constraints constraints = [A_wall >= 2*h*w + 2*h*d, A_floor >= w*d, h/w >= alpha, h/w <= beta, d/w >= gamma, d/w <= delta] #Objective function V = h*w*d objective = 1/V #To maximize V, we minimize its reciprocal # Formulate the GP gp = GP(objective, constraints) # Solve the GP sol = gp.solve() # Print results table print sol.table()Water Tank¶
from gpkit import Variable, VectorVariable, GP M = Variable("M", 100, "kg", "Mass of Water in the Tank") rho = Variable("\\rho", 1000, "kg/m^3", "Density of Water in the Tank") A = Variable("A", "m^2", "Surface Area of the Tank") V = Variable("V", "m^3", "Volume of the Tank") d = VectorVariable(3, "d", "m", "Dimension Vector") constraints = (A >= 2*(d[0]*d[1] + d[0]*d[2] + d[1]*d[2]), V == d[0]*d[1]*d[2], M == V*rho ) gp = GP(A, constraints) sol = gp.solve(printing=False) print sol(A) print sol(V) print sol(d)iPython Notebook Examples¶
Also available on nbviewer.
BOX¶
Maximize volume while limiting the aspect ratios and area of the sides.
This is a simple demonstration of gpkit, based on an example from A Tutorial on Geometric Programming by Boyd et al..
Set up the modelling environment¶
First we’ll to import GPkit and turn on \(\LaTeX\) printing for GPkit variables and equations.
import gpkit import gpkit.interactive gpkit.interactive.init_printing()Now we declare the optimisation parameters:
alpha = gpkit.Variable("\\alpha", 2, "-", "lower limit, wall aspect ratio") beta = gpkit.Variable("\\beta", 10, "-", "upper limit, wall aspect ratio") gamma = gpkit.Variable("\\gamma", 2, "-", "lower limit, floor aspect ratio") delta = gpkit.Variable("\\delta", 10, "-", "upper limit, floor aspect ratio") A_wall = gpkit.Variable("A_{wall}", 200, "m^2", "upper limit, wall area") A_floor = gpkit.Variable("A_{floor}", 50, "m^2", "upper limit, floor area")Next, we declare the decision variables:
var = gpkit.Variable # a convenient shorthand h = var("h", "m", "height") w = var("w", "m", "width") d = var("d", "m", "depth")Then we form equations of the system:
V = h*w*d constraints = [A_wall >= 2*h*w + 2*h*d, A_floor >= w*d, h/w >= alpha, h/w <= beta, d/w >= gamma, d/w <= delta]Formulate the optimisation problem¶
Here we write the optimisation problem as a standard form geometric program. Note that by putting \(\tfrac{1}{V}\) as our cost function, we are maximizing \(V\).
gp = gpkit.GP(1/V, constraints)Now we can check that our equations are correct by using the built-in latex printing.
gp\[\begin{split}\begin{array}[ll] \text{} \text{minimize} & \frac{1}{d h w}\mathrm{\left[ \tfrac{1}{m^{3}} \right]} \\ \text{subject to} & A_{wall} \geq 2d h + 2h w \\ & A_{floor} \geq d w \\ & \alpha \leq \frac{h}{w} \\ & \beta \geq \frac{h}{w} \\ & \gamma \leq \frac{d}{w} \\ & \delta \geq \frac{d}{w} \\ \text{substituting} & A_{floor} = 50 \\ & A_{wall} = 200 \\ & \alpha = 2 \\ & \beta = 10 \\ & \delta = 10 \\ & \gamma = 2 \\ \end{array}\end{split}\]That looks fine, but let’s change \(A_{floor}\) to \(100\), just for fun.
Note that replace=True will be needed, since \(A_{floor}\) has already been substituted in.
gp.sub(A_floor, 100, replace=True) # (var, val) substitution syntaxAnd heck, why not change \(\alpha\) and \(\gamma\) to \(1\) while we’re at it?
gp.sub({alpha: 1, gamma: 1}, replace=True) # ({var1: val1, var2: val2}) substitution syntaxNow check that those changes took:
gp\[\begin{split}\begin{array}[ll] \text{} \text{minimize} & \frac{1}{d h w}\mathrm{\left[ \tfrac{1}{m^{3}} \right]} \\ \text{subject to} & A_{wall} \geq 2d h + 2h w \\ & A_{floor} \geq d w \\ & \alpha \leq \frac{h}{w} \\ & \beta \geq \frac{h}{w} \\ & \gamma \leq \frac{d}{w} \\ & \delta \geq \frac{d}{w} \\ \text{substituting} & A_{floor} = 100 \\ & A_{wall} = 200 \\ & \alpha = 1 \\ & \beta = 10 \\ & \delta = 10 \\ & \gamma = 1 \\ \end{array}\end{split}\]Looks good!
Analyse results¶
print sol.table()0.0025981 : Cost | Free variables d : 11.5 [m] depth h : 5.77 [m] height w : 5.77 [m] width | | Constants A_{floor} : 100 [m**2] upper limit, floor area A_{wall} : 200 [m**2] upper limit, wall area alpha : 1 [-] lower limit, wall aspect ratio beta : 10 [-] upper limit, wall aspect ratio delta : 10 [-] upper limit, floor aspect ratio gamma : 1 [-] lower limit, floor aspect ratio | | Constant sensitivities A_{wall} : -1.5 [-] upper limit, wall area alpha : 0.5 [-] lower limit, wall aspect ratio |Hmm, why didn’t \(A_{floor}\) show up in the sensitivities list?
sol.senssubinto(A_floor)-3.4698494246624108e-09Its sensitivity is tiny; changing it near this value doesn’t affect the cost at all; that constraint is loose! Let’s sweep over a range of \(A_{floor}\) values to figure out where it becomes loose.
Sweep and plot results¶
Import the plotting library matplotlib and the math library numpy:
%matplotlib inline %config InlineBackend.figure_format = 'retina' # for high-DPI displays import matplotlib.pyplot as plt import numpy as npSolve for values of \(A_{floor}\) from 10 to 100 using a “sweep” substitution:
gp.sub(A_floor, ('sweep', np.linspace(10, 100, 50)), replace=True) sol = gp.solve() print sol.table()Using solver 'cvxopt' Sweeping 1 variables over 50 passes Solving took 0.547 seconds 0.0032211 : Cost (average of 50 values) | Free variables (average) d : 8.33 [m] depth h : 7.74 [m] height w : 5.69 [m] width | | Constants (average) A_{floor} : 55 [m**2] upper limit, floor area A_{wall} : 200 [m**2] upper limit, wall area alpha : 1 [-] lower limit, wall aspect ratio beta : 10 [-] upper limit, wall aspect ratio delta : 10 [-] upper limit, floor aspect ratio gamma : 1 [-] lower limit, floor aspect ratio | | Constant sensitivities (average) A_{floor} : -0.274 [-] upper limit, floor area A_{wall} : -1.23 [-] upper limit, wall area alpha : 0.226 [-] lower limit, wall aspect ratio |It seems we got some sensitivity out of \(A_{floor}\) on average over these points; let’s plot it:
plt.plot(sol(A_floor), sol(d), linewidth=1, alpha=0.5) plt.plot(sol(A_floor), sol(h), linewidth=1, alpha=0.5) plt.plot(sol(A_floor), sol(w), '--', linewidth=2, alpha=0.5) plt.legend(['depth', 'height', 'width']) plt.ylabel("Optimal dimensions [m]") _ = plt.xlabel("A_floor [m^2]") # the _ catches the returned label object, since we don't need it
There’s an interesting elbow when \(A_{floor} \approx 50 \mathrm{\ m^2}\).
Interactive analysis¶
Let’s investigate it with the cadtoons library. Running cadtoon.py box.svg in this folder creates an interactive SVG graphic for us.
First, import the functions to display HTML in iPython Notebook, and the ractivejs library.
from IPython.display import HTML, display from string import TemplateThen write controls that link the optimization to the animation. It’s generally very helpful to play around with with writing constraints in the cadtoons sandbox before copying the javascript code from there to iPython.
ractor = Template(""" var w = $w/10, h = $h/10, d = $d/10""") def ractorfn(sol): return ractor.substitute(w=sol(w), d=sol(d), h=sol(h)) constraints=""" var dw = 50*(w-1), dh = 50*(h-1), dd = 50*(d-1) r.plan.scalex = w r.plan.scaley = d r.top.scalex = w r.top.scaley = h r.top.y = -dh -dd r.bottom.scalex = w r.bottom.scaley = h r.bottom.y = dh + dd r.left.scalex = h r.left.scaley = d r.left.x = -dh - dw r.right.scalex = h r.right.scaley = d r.right.x = dh + dw """ def ractivefn(gp): sol = gp.solution live = "<script>" + ractorfn(sol) + constraints + "</script>" display(HTML(live)) # if you enable the line below, you can try navigating the sliders by sensitivities # print sol.table(["cost", "sensitivities"])Now that the display window has been created, let’s use an iPython widget to explore the “design space” of our box. This widget runs the gpkit solver every time a slider is changed, so you only solve the points you see.
with open("box.gpkit", 'r') as file: display(HTML(file.read())) display(HTML(display(HTML("<style> #ractivecontainer" "{position:absolute; height: 0;" "right: 0; top: 5em;} </style>"))))<IPython.core.display.HTML at 0x7fd2dce97850>gpkit.interactive.widget(gp, ractivefn, { "A_{floor}": (10, 1000, 10), "A_{wall}": (10, 1000, 10), "\\delta": (0.1, 20, 0.1), "\\gamma": (0.1, 20, 0.1), "\\alpha": (0.1, 20, 0.1), "\\beta": (0.1, 20, 0.1), })Now that we’ve found a good range to explore, we’ll make and interactive javascript widget that presolves all options and works in a iPython notebook hosted by nbviewer.
gpkit.interactive.jswidget(gp, ractorfn, constraints, { "A_{floor}": (10, 100, 10), "A_{wall}": (10, 210, 20), "\\delta": (1, 10, 3), "\\gamma": (1, 10, 3), "\\alpha": (1, 10, 3), "\\beta": (1, 10, 3), })This concludes the Box example. Try playing around with the sliders up above until you’re bored of this simple system; then check out one of the other examples. Thanks for reading!
Import CSS for nbviewer¶
If you have a local iPython stylesheet installed, this will add it to the iPython Notebook:
from IPython import utils from IPython.core.display import HTML import os def css_styling(): """Load default custom.css file from ipython profile""" base = utils.path.get_ipython_dir() styles = ("<style>\n%s\n</style>" % open(os.path.join(base,'profile_default/static/custom/custom.css'),'r').read()) return HTML(styles) css_styling()
AIRPLANE FUEL¶
Minimize fuel needed for a plane that can sprint and land quickly.
Set up the modelling environment¶
First we’ll to import GPkit and turn on \(\LaTeX\) printing for GPkit variables and equations.
import numpy as np import gpkit import gpkit.interactive gpkit.interactive.init_printing()declare constants¶
mon = gpkit.Variable vec = gpkit.VectorVariable N_lift = mon("N_{lift}", 6.0, "-", "Wing loading multiplier") pi = mon("\\pi", np.pi, "-", "Half of the circle constant") sigma_max = mon("\\sigma_{max}", 250e6, "Pa", "Allowable stress, 6061-T6") sigma_maxshear = mon("\\sigma_{max,shear}", 167e6, "Pa", "Allowable shear stress") g = mon("g", 9.8, "m/s^2", "Gravitational constant") w = mon("w", 0.5, "-", "Wing-box width/chord") r_h = mon("r_h", 0.75, "-", "Wing strut taper parameter") f_wadd = mon("f_{wadd}", 2, "-", "Wing added weight fraction") W_fixed = mon("W_{fixed}", 14.7e3, "N", "Fixed weight") C_Lmax = mon("C_{L,max}", 1.5, "-", "Maximum C_L, flaps down") rho = mon("\\rho", 0.91, "kg/m^3", "Air density, 3000m") rho_sl = mon("\\rho_{sl}", 1.23, "kg/m^3", "Air density, sea level") rho_alum = mon("\\rho_{alum}", 2700, "kg/m^3", "Density of aluminum") mu = mon("\\mu", 1.69e-5, "kg/m/s", "Dynamic viscosity, 3000m") e = mon("e", 0.95, "-", "Wing spanwise efficiency") A_prop = mon("A_{prop}", 0.785, "m^2", "Propeller disk area") eta_eng = mon("\\eta_{eng}", 0.35, "-", "Engine efficiency") eta_v = mon("\\eta_v", 0.85, "-", "Propeller viscous efficiency") h_fuel = mon("h_{fuel}", 42e6, "J/kg", "fuel heating value") V_sprint_reqt = mon("V_{sprintreqt}", 150, "m/s", "sprint speed requirement") W_pay = mon("W_{pay}", 500*9.81, "N") R_min = mon("R_{min}", 1e6, "m", "Minimum airplane range") V_stallmax = mon("V_{stall,max}", 40, "m/s", "Stall speed") # sweep variables R_min = mon("R_{min}", np.linspace(1e6, 1e7, 10), "m", "Minimum airplane range") V_stallmax = mon("V_{stall,max}", np.linspace(30, 50, 10), "m/s", "Stall speed")declare free variables¶
V = vec(3, "V", "m/s", "Flight speed") C_L = vec(3, "C_L", "-", "Wing lift coefficent") C_D = vec(3, "C_D", "-", "Wing drag coefficent") C_Dfuse = vec(3, "C_{D_{fuse}}", "-", "Fuselage drag coefficent") C_Dp = vec(3, "C_{D_p}", "-", "drag model parameter") C_Di = vec(3, "C_{D_i}", "-", "drag model parameter") T = vec(3, "T", "N", "Thrust force") Re = vec(3, "Re", "-", "Reynold's number") W = vec(3, "W", "N", "Aircraft weight") eta_i = vec(3, "\\eta_i", "-", "Aircraft efficiency") eta_prop = vec(3, "\\eta_{prop}", "-") eta_0 = vec(3, "\\eta_0", "-") W_fuel = vec(2, "W_{fuel}", "N", "Fuel weight") z_bre = vec(2, "z_{bre}", "-") S = mon("S", "m^2", "Wing area") R = mon("R", "m", "Airplane range") A = mon("A", "-", "Aspect Ratio") I_cap = mon("I_{cap}", "m^4", "Spar cap area moment of inertia per unit chord") M_rbar = mon("\\bar{M}_r", "-") P_max = mon("P_{max}", "W") V_stall = mon("V_{stall}", "m/s") nu = mon("\\nu", "-") p = mon("p", "-") q = mon("q", "-") tau = mon("\\tau", "-") t_cap = mon("t_{cap}", "-") t_web = mon("t_{web}", "-") W_cap = mon("W_{cap}", "N") W_zfw = mon("W_{zfw}", "N", "Zero fuel weight") W_eng = mon("W_{eng}", "N") W_mto = mon("W_{mto}", "N", "Maximum takeoff weight") W_pay = mon("W_{pay}", "N") W_tw = mon("W_{tw}", "N") W_web = mon("W_{web}", "N") W_wing = mon("W_{wing}", "N")Let’s check that the vector constraints are working:
W == 0.5*rho*C_L*S*V**2\[\begin{split}\begin{bmatrix}{W}_{0} = 0.5S \rho {V}_{0}^{2} {C_L}_{0} & {W}_{1} = 0.5S \rho {V}_{1}^{2} {C_L}_{1} & {W}_{2} = 0.5S \rho {V}_{2}^{2} {C_L}_{2}\end{bmatrix}\end{split}\]Looks good!
Form the optimization problem¶
In the 3-element vector variables, indexs 0, 1, and 2 are the outbound, return and sprint flights.
steady_level_flight = (W == 0.5*rho*C_L*S*V**2, T >= 0.5*rho*C_D*S*V**2, Re == (rho/mu)*V*(S/A)**0.5) landing_fc = (W_mto <= 0.5*rho_sl*V_stall**2*C_Lmax*S, V_stall <= V_stallmax) sprint_fc = (P_max >= T[2]*V[2]/eta_0[2], V[2] >= V_sprint_reqt) drag_model = (C_D >= (0.05/S)*gpkit.units.m**2 +C_Dp + C_L**2/(pi*e*A), 1 >= (2.56*C_L**5.88/(Re**1.54*tau**3.32*C_Dp**2.62) + 3.8e-9*tau**6.23/(C_L**0.92*Re**1.38*C_Dp**9.57) + 2.2e-3*Re**0.14*tau**0.033/(C_L**0.01*C_Dp**0.73) + 1.19e4*C_L**9.78*tau**1.76/(Re*C_Dp**0.91) + 6.14e-6*C_L**6.53/(Re**0.99*tau**0.52*C_Dp**5.19))) propulsive_efficiency = (eta_0 <= eta_eng*eta_prop, eta_prop <= eta_i*eta_v, 4*eta_i + T*eta_i**2/(0.5*rho*V**2*A_prop) <= 4) # 4th order taylor approximation for e^x z_bre_sum = 0 for i in range(1,5): z_bre_sum += z_bre**i/np.math.factorial(i) range_constraints = (R >= R_min, z_bre >= g*R*T[:2]/(h_fuel*eta_0[:2]*W[:2]), W_fuel/W[:2] >= z_bre_sum) weight_relations = (W_pay >= 500*g*gpkit.units.kg, W_tw >= W_fixed + W_pay + W_eng, W_zfw >= W_tw + W_wing, W_eng >= 0.0372*P_max**0.8083 * gpkit.units.parse_expression('N/W^0.8083'), W_wing/f_wadd >= W_cap + W_web, W[0] >= W_zfw + W_fuel[1], W[1] >= W_zfw, W_mto >= W[0] + W_fuel[0], W[2] == W[0]) wing_structural_model = (2*q >= 1 + p, p >= 1.9, tau <= 0.15, M_rbar >= W_tw*A*p/(24*gpkit.units.N), .92**2/2.*w*tau**2*t_cap >= I_cap * gpkit.units.m**-4 + .92*w*tau*t_cap**2, 8 >= N_lift*M_rbar*A*q**2*tau/S/I_cap/sigma_max * gpkit.units.parse_expression('Pa*m**6'), 12 >= A*W_tw*N_lift*q**2/tau/S/t_web/sigma_maxshear, nu**3.94 >= .86*p**-2.38 + .14*p**0.56, W_cap >= 8*rho_alum*g*w*t_cap*S**1.5*nu/3/A**.5, W_web >= 8*rho_alum*g*r_h*tau*t_web*S**1.5*nu/3/A**.5 )eqns = (weight_relations + range_constraints + propulsive_efficiency + drag_model + steady_level_flight + landing_fc + sprint_fc + wing_structural_model) gp = gpkit.GP(W_fuel.sum(), eqns)Design a hundred airplanes¶
sol = gp.solve()Using solver 'cvxopt' Sweeping 2 variables over 100 passes Solving took 9.96 secondsThe “local model” is the power-law tangent to the Pareto frontier, gleaned from sensitivities.
sol["local_model"][0].sub(gp.substitutions)\[0.006484\frac{R_{min}}{V_{stall,max}^{0.59}}\]plot design frontiers¶
%matplotlib inline %config InlineBackend.figure_format = 'retina' plot_frontiers = gpkit.interactive.plot_frontiers plot_frontiers(gp, [V[0], V[1], V[2]]) plot_frontiers(gp, [S, W_zfw, P_max]) plot_frontiers(gp, ['S{\\rho_{alum}}', 'S{h_{fuel}}', 'S{A_{prop}}'])
Interactive analysis¶
Let’s investigate it with the cadtoons library. Running cadtoon.py flightconditions.svg in this folder creates an interactive SVG graphic for us.
First, import the functions to display HTML in iPython Notebook, and the ractivejs library.
from IPython.display import HTML, display from string import Templateractor = Template(""" var W_eng = $W_eng, lam = $lam r.shearinner.scalex = 1-$tcap*10 r.shearinner.scaley = 1-$tweb*100 r.airfoil.scaley = $tau/0.13 r.fuse.scalex = $W_fus/24000 r.wing.scalex = $b/2/14 r.wing.scaley = $cr*1.21 """) def ractorfn(sol): return ractor.substitute(lam = 0.5*(sol(p) - 1), b = sol((S*A)**0.5), cr = sol(2/(1+0.5*(sol(p) - 1))*(S/A)**0.5), tcap = sol(t_cap/tau), tweb = sol(t_web/w), tau = sol(tau), W_eng = sol(W_eng), W_fus = sol(W_mto) - sol(W_wing) - sol(W_eng)) constraints = """ r.engine1.scale = Math.pow(W_eng/3000, 2/3) r.engine2.scale = Math.pow(W_eng/3000, 2/3) r.engine1.y = 6*lam r.engine2.y = 6*lam r.wingrect.scaley = 1-lam r.wingrect.y = -6 + 5*lam r.wingtaper.scaley = lam r.wingtaper.y = 5*lam """ def ractivefn(gp): sol = gp.solution live = "<script>" + ractorfn(sol) + constraints + "</script>" display(HTML(live)) # if you enable the line below, you can try navigating the sliders by sensitivities # print sol.table(["cost", "sensitivities"])with open("flightconditions.gpkit", 'r') as file: display(HTML(file.read())) display(HTML("<style> #ractivecontainer" "{position:absolute; height: 0;" "right: 0; top: -6em;} </style>"))gpkit.interactive.widget(gp, ractivefn, {"V_{stall,max}": (20, 50, 1), "R_{min}": (1e6, 1e7, 0.5e6)})gpkit.interactive.jswidget(gp, ractorfn, constraints, {"V_{stall,max}": (20, 50, 3), "R_{min}": (1e6, 1e7, 1e6)})This concludes the Box example. Try playing around with the sliders up above until you’re bored of this system; then check out one of the other examples. Thanks for reading!
Import CSS for nbviewer¶
If you have a local iPython stylesheet installed, this will add it to the iPython Notebook:
from IPython import utils from IPython.core.display import HTML import os def css_styling(): """Load default custom.css file from ipython profile""" base = utils.path.get_ipython_dir() styles = "<style>\n%s\n</style>" % (open(os.path.join(base,'profile_default/static/custom/custom.css'),'r').read()) return HTML(styles) css_styling()Glossary¶
For an alphabetical listing of all commands, check out the Index
The GPkit Package¶
Lightweight GP Modeling Package
For examples please see the examples folder.
Requirements¶
numpy MOSEK or CVXOPT scipy(optional): for complete sparse matrix support sympy(optional): for latex printing in iPython Notebook
Attributes¶
- settings : dict
- Contains settings loaded from ./env/settings
- gpkit.disableUnits()¶
Disables units support in a particular instance of GPkit.
Posynomials created after this is run are incompatible with those created before. If gpkit is imported multiple times, this needs to be run each time.
- gpkit.enableUnits()¶
Enables units support in a particular instance of GPkit.
Posynomials created after this is run are incompatible with those created before. If gpkit is imported multiple times, this needs to be run each time.
Subpackages¶
Submodules¶
gpkit.geometric_program¶
Module for creating GP instances.
Example¶
>>> gp = gpkit.GP(cost, constraints, substitutions)
- class gpkit.geometric_program.GP(cost, constraints, substitutions={}, solver=None, options={})¶
Bases: gpkit.model.Model
Holds a model and cost function for passing to solvers.
- cost : Constraint
- Posynomial to minimize when solving
- constraints : list of (lists of) Constraints
- Constraints to maintain when solving (MonoEQConstraints will be turned into <= and >= constraints)
- substitutions : dict {varname: float or int} (optional)
- Substitutions to be applied before solving (including sweeps)
- solver : str (optional)
- Name of solver to use
- options : dict (optional)
- Options to pass to solver
>>> gp = gpkit.GP( # minimize 0.5*rho*S*C_D*V**2, [ # subject to Re <= (rho/mu)*V*(S/A)**0.5, C_f >= 0.074/Re**0.2, W <= 0.5*rho*S*C_L*V**2, W <= 0.5*rho*S*C_Lmax*V_min**2, W >= W_0 + W_w, W_w >= W_w_surf + W_w_strc, C_D >= C_D_fuse + C_D_wpar + C_D_ind ], substitutions) >>> gp.solve()
- solve(solver=None, printing=True, skipfailures=False)¶
Solves a GP and returns the solution.
- printing : bool (optional)
- If True (default), then prints out solver used and time to solve.
- solution : dict
- A dictionary containing the optimal values for each free variable.
- class gpkit.geometric_program.GPSolutionArray¶
Bases: gpkit.small_classes.DictOfLists
DictofLists extended with posynomial substitution.
- senssubinto(p)¶
Returns array of each solution’s sensitivity substituted into p
Returns only scalar values.
- subinto(p)¶
Returns numpy array of each solution substituted into p.
- table(tables=['cost', 'free_variables', 'constants', 'sensitivities'])¶
gpkit.model¶
Module for creating models.
Currently these are only used for GP instances, but they may be further abstractable.
gpkit.nomials¶
- class gpkit.nomials.Constraint(p1, p2)¶
Bases: gpkit.nomials.Posynomial
TODO: Add docstring
- class gpkit.nomials.MonoEQConstraint(p1, p2)¶
Bases: gpkit.nomials.Constraint
TODO: Add docstring
- class gpkit.nomials.Monomial(exps=None, cs=1, var_locs=None, allow_negative=False, **descr)¶
Bases: gpkit.nomials.Posynomial
TODO: Add docstring
- class gpkit.nomials.Posynomial(exps=None, cs=1, var_locs=None, allow_negative=False, **descr)¶
Bases: object
A representation of a posynomial.
- exps: tuple of dicts
- Exponent dicts for each monomial term
- cs: tuple
- Coefficient values for each monomial term
- var_locs: dict
- mapping from variable name to list of indices of monomial terms that variable appears in
Posynomial (if the input has multiple terms) Monomial (if the input has one term)
- descr(descr)¶
- sub(substitutions, val=None, allow_negative=False)¶
- subcmag(substitutions, val=None)¶
- to(arg)¶
- class gpkit.nomials.VarKey(k=None, **kwargs)¶
Bases: object
Used in the declaration of Posynomials to correspond to each ‘variable name’.
- k : object (usually str)
- The variable’s name attribute is derived from str(k).
- **kwargs :
- Any additional attributes, which become the descr attribute (a dict).
VarKey with the given name and descr.
- new_unnamed_id = <method-wrapper 'next' of itertools.count object at 0x7fcf836db758>¶
- class gpkit.nomials.Variable(*args, **descr)¶
Bases: gpkit.nomials.Monomial
- class gpkit.nomials.VectorVariable¶
Bases: gpkit.posyarray.PosyArray
A described vector of singlet Monomials.
- length : int
- Length of vector.
- *args : list
- may contain “name” (Strings)
and/or “label” (Strings)
- **descr : dict
- VarKey description
PosyArray of Monomials, each containing a VarKey with name ‘$name_{i}’, where $name is the vector’s name and i is the VarKey’s index.
gpkit.posyarray¶
Module for creating PosyArray instances.
Example¶
>>> x = gpkit.Monomial('x') >>> px = gpkit.PosyArray([1, x, x**2])
- class gpkit.posyarray.PosyArray¶
Bases: numpy.ndarray
A Numpy array with elementwise inequalities and substitutions.
input_array : array-like
>>> px = gpkit.PosyArray([1, x, x**2])
- c¶
- outer(other)¶
Returns the array and argument’s outer product.
- prod()¶
Returns the product of the array.
- sub(subs, val=None, allow_negative=False)¶
Substitutes into the array
- sum()¶
Returns the sum of the array.
gpkit.small_classes¶
- class gpkit.small_classes.CootMatrix¶
Bases: gpkit.small_classes.CootMatrix
A very simple sparse matrix representation.
- append(i, j, x)¶
- shape = (None, None)¶
- tocoo()¶
Converts to a Scipy sparse coo_matrix
- tocsc()¶
- tocsr()¶
- todense()¶
- todia()¶
- todok()¶
- update_shape()¶
- gpkit.small_classes.CootMatrixTuple¶
alias of CootMatrix
- class gpkit.small_classes.DictOfLists¶
Bases: dict
A hierarchy of dicionaries, with lists at the bottom.
- append(sol)¶
Appends a dict (of dicts) of lists to all held lists.
- atindex(i)¶
Indexes into each list independently.
- toarray(shape=None)¶
Converts all lists into arrays.
- class gpkit.small_classes.HashVector¶
Bases: dict
A simple, sparse, string-indexed immutable vector. Inherits from dict.
The HashVector class supports element-wise arithmetic: any undeclared variables are assumed to have a value of zero.
arg : iterable
>>> x = gpkit.nomials.Monomial('x') >>> exp = gpkit.small_classes.HashVector({x: 2})
- class gpkit.small_classes.PosyTuple¶
Bases: tuple
PosyTuple(exps, cs, var_locs, substitutions)
- cs¶
Alias for field number 1
- exps¶
Alias for field number 0
- substitutions¶
Alias for field number 3
- var_locs¶
Alias for field number 2
- gpkit.small_classes.append_dict(i, o)¶
Recursviely travels dict o and appends items found in i.
- gpkit.small_classes.enlist_dict(i, o)¶
Recursviely copies dict i into o, placing non-dict items into lists.
- gpkit.small_classes.enray_dict(i, o)¶
Recursively turns lists into numpy arrays.
- gpkit.small_classes.index_dict(idx, i, o)¶
Recursviely travels dict i, placing items at idx into dict o.
gpkit.small_scripts¶
- gpkit.small_scripts.flatten(ible, classes)¶
Flatten an iterable that contains other iterables
- l : Iterable
- Top-level container
- out : list
- List of all objects found in the nested iterables
- TypeError
- If an object is found whose class was not in classes
- gpkit.small_scripts.invalid_types_for_oper(oper, a, b)¶
Raises TypeError for unsupported operations.
- gpkit.small_scripts.is_sweepvar(sub)¶
Determines if a given substitution indicates a sweep.
- gpkit.small_scripts.latex_num(c)¶
- gpkit.small_scripts.locate_vars(exps)¶
From exponents form a dictionary of which monomials each variable is in.
- gpkit.small_scripts.mag(c)¶
Return magnitude of a Number or Quantity
- gpkit.small_scripts.results_table(data, title, senss=False)¶
- gpkit.small_scripts.sort_and_simplify(exps, cs)¶
Reduces the number of monomials, and casts them to a sorted form.
- gpkit.small_scripts.unitstr(v, into='%s', options='~')¶
gpkit.substitution¶
Module containing the substitution function
- gpkit.substitution.substitution(var_locs, exps, cs, substitutions, val=None)¶
Efficient substituton into a list of monomials.
- var_locs : dict
- Dictionary of monomial indexes for each variable.
- exps : dict
- Dictionary of variable exponents for each monomial.
- cs : list
- Coefficients each monomial.
- substitutions : dict
- Substitutions to apply to the above.
- val : number (optional)
- Used to substitute singlet variables.
- gpkit.substitution.vectorsub(subs, var, sub, varset)¶
Vectorized substitution via vecmons and Variables.
Citing GPkit¶
If you use GPkit, please cite it using the following bibtex:
@Misc{gpkit, author = {MIT Convex Optimization Group}, title = {GPkit}, howpublished = {\url{http://github.com/convexopt/gpkit}}, year = {2015}, note = {Version 0.1} }