Advanced Commands

Choice Variables

If MOSEK 9 is installed, GPkit supports discretized free variables with the mosek_conif solver. Choice variables are free in the sense of having multiple possible choices, but discretized in the sense of having a limited set of possible solutions.

"Example choice variable usage"
import numpy as np
from gpkit import Variable, Model

x = Variable("x", choices=range(1, 4))
num = Variable("numerator", np.linspace(0.5, 7, 11))

m = Model(x + num/x)
sol = m.solve(verbosity=0)

print(sol.table())

If solved with the mosek_conif solver, the script above will print:

Optimal Cost
------------
 [ 1.5       2.15      2.8       3.22     ... ]

~~~~~~~~
WARNINGS
~~~~~~~~
No Dual Solution
----------------
This model has the discretized choice variables [x] and hence no dual
solution. You can fix those variables to their optimal value and get
sensitivities to the resulting continuous problem by updating your model's
substitions with `sol["choicevariables"]`.
~~~~~~~~

Swept Variables
---------------
numerator : [ 0.5
              1.15
              1.8
              2.45
              3.1
              3.75
              4.4
              5.05
              5.7
              6.35
              7         ]

Free Variables
--------------
x : [ 1
      1
      1
      2
      2
      2
      2
      2
      2
      3
      3         ]

Note that the optimal values for x are discretized, clicking from 1 to 2 to 3 during the sweep, and that the solution has no dual variables.

Derived Variables

Evaluated Fixed Variables

Some fixed variables may be derived from the values of other fixed variables. For example, air density, viscosity, and temperature are functions of altitude. These can be represented by a substitution or value that is a one-argument function accepting model.substitutions (for details, see Substitutions below).

"Example pre-solve evaluated fixed variable"
from gpkit import Variable, Model, units

# code from t_GPSubs.test_calcconst in tests/t_sub.py
x = Variable("x", "hours")
t_day = Variable("t_{day}", 12, "hours")
t_night = Variable("t_{night}",
                   lambda c: 1*units.day - c(t_day), "hours")

# note that t_night has a function as its value
m = Model(x, [x >= t_day, x >= t_night])
sol = m.solve(verbosity=0)
assert sol["variables"][t_night] == 12

# call substitutions
m.substitutions.update({t_day: ("sweep", [8, 12, 16])})
sol = m.solve(verbosity=0)
assert (sol["variables"][t_night] == [16, 12, 8]).all()

These functions are automatically differentiated with the ad package to provide more accurate sensitivities. In some cases may require using functions from the ad.admath instead of their python or numpy equivalents; the ad documentation contains details on how to do this.

Evaluated Free Variables

Some free variables may be evaluated from the values of other (non-evaluated) free variables after the optimization is performed. For example, if the efficiency \(\nu\) of a motor is not a GP-compatible variable, but \((1-\nu)\) is a valid GP variable, then \(\nu\) can be calculated after solving. These evaluated free variables can be represented by a Variable with evalfn metadata. When constructing an evalfn, remember that square-bracket access to variables pulls out magnitudes: use round-bracket access (i.e. v(var)) to ensure unit correctness.

"Example post-solve evaluated variable"
from gpkit import Variable, Model

# code from t_constraints.test_evalfn in tests/t_sub.py
x = Variable("x")
x2 = Variable("x^2", evalfn=lambda v: v(x)**2)
m = Model(x, [x >= 2])
m.unique_varkeys = set([x2.key])
sol = m.solve(verbosity=0)
assert abs(sol(x2) - 4) <= 1e-4

Note that this variable should not be used in constructing your model! For evaluated variables that can be used during a solution, see Sequential Geometric Programs.

Sweeps

Sweeps are useful for analyzing tradeoff surfaces. A sweep “value” is an Iterable of numbers, e.g. [1, 2, 3]. The simplest way to sweep a model is to call model.sweep({sweepvar: sweepvalues}), which will return a solution array but not change the model’s substitutions dictionary. If multiple sweepvars are given, the method will run them all as independent one-dimensional sweeps and return a list of one solution per sweep. The method model.autosweep({sweepvar: (start, end)}, tol=0.01) behaves very similarly, except that only the bounds of the sweep need be specified and the region in betwen will be swept to a maximum possible error of tol in the log of the cost. For details see 1D Autosweeps below.

Sweep Substitutions

Alternatively, or to sweep a higher-dimensional grid, Variables can swept with a substitution value takes the form ('sweep', Iterable), such as ('sweep', np.linspace(1e6, 1e7, 100)). During variable declaration, giving an Iterable value for a Variable is assumed to be giving it a sweep value: for example, x = Variable("x", [1, 2, 3]) will sweep x over three values.

Vector variables may also be substituted for: {y: ("sweep" ,[[1, 2], [1, 2], [1, 2]])} will sweep \(y\ \forall~y_i\in\left\{1,2\right\}\). These sweeps cannot be specified during Variable creation.

A Model with sweep substitutions 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.

1D Autosweeps

If you’re only sweeping over a single variable, autosweeping lets you specify a tolerance for cost error instead of a number of exact positions to solve at. GPkit will then search the sweep segment for a locally optimal number of sweeps that can guarantee a max absolute error on the log of the cost.

Accessing variable and cost values from an autosweep is slightly different, as can be seen in this example:

"Show autosweep_1d functionality"
import pickle
import numpy as np
import gpkit
from gpkit import units, Variable, Model
from gpkit.tools.autosweep import autosweep_1d
from gpkit.small_scripts import mag

A = Variable("A", "m**2")
l = Variable("l", "m")

m1 = Model(A**2, [A >= l**2 + units.m**2])
tol1 = 1e-3
bst1 = autosweep_1d(m1, tol1, l, [1, 10], verbosity=0)
print("Solved after %2i passes, cost logtol +/-%.3g" % (bst1.nsols, bst1.tol))
# autosweep solution accessing
l_vals = np.linspace(1, 10, 10)
sol1 = bst1.sample_at(l_vals)
print("values of l: %s" % l_vals)
print("values of A: [%s] %s" %
      (" ".join("% .1f" % n for n in sol1("A").magnitude), sol1("A").units))
cost_estimate = sol1["cost"]
cost_lb, cost_ub = sol1.cost_lb(), sol1.cost_ub()
print("cost lower bound:\n%s\n" % cost_lb)
print("cost estimate:\n%s\n" % cost_estimate)
print("cost upper bound:\n%s\n" % cost_ub)
# you can evaluate arbitrary posynomials
np.testing.assert_allclose(mag(2*sol1(A)), mag(sol1(2*A)))
assert (sol1["cost"] == sol1(A**2)).all()
# the cost estimate is the logspace mean of its upper and lower bounds
np.testing.assert_allclose((np.log(mag(cost_lb)) + np.log(mag(cost_ub)))/2,
                           np.log(mag(cost_estimate)))
# save autosweep to a file and retrieve it
bst1.save("autosweep.pkl")
bst1_loaded = pickle.load(open("autosweep.pkl", "rb"))

# this problem is two intersecting lines in logspace
m2 = Model(A**2, [A >= (l/3)**2, A >= (l/3)**0.5 * units.m**1.5])
tol2 = {"mosek_cli": 1e-6, "mosek_conif": 1e-6,
        "cvxopt": 1e-7}[gpkit.settings["default_solver"]]
# test Model method
sol2 = m2.autosweep({l: [1, 10]}, tol2, verbosity=0)
bst2 = sol2.bst
print("Solved after %2i passes, cost logtol +/-%.3g" % (bst2.nsols, bst2.tol))
print("Table of solutions used in the autosweep:")
print(bst2.solarray.table())

If you need access to the raw solutions arrays, the smallest simplex tree containing any given point can be gotten with min_bst = bst.min_bst(val), the extents of that tree with bst.bounds and solutions of that tree with bst.sols. More information is in help(bst).

Tight ConstraintSets

Tight ConstraintSets will warn if any inequalities they contain are not tight (that is, the right side does not equal the left side) after solving. This is useful when you know that a constraint should be tight for a given model, but representing it as an equality would be non-convex.

"Example Tight ConstraintSet usage"
from gpkit import Variable, Model
from gpkit.constraints.tight import Tight

Tight.reltol = 1e-2  # set the global tolerance of Tight
x = Variable('x')
x_min = Variable('x_{min}', 2)
m = Model(x, [Tight([x >= 1], reltol=1e-3),  # set the specific tolerance
              x >= x_min])
m.solve(verbosity=0)  # prints warning

Loose ConstraintSets

Loose ConstraintSets will warn if any GP-compatible constraints they contain are not loose (that is, their sensitivity is above some threshold after solving). This is useful when you want a constraint to be inactive for a given model because it represents an important model assumption (such as a fit only valid over a particular interval).

"Example Loose ConstraintSet usage"
from gpkit import Variable, Model
from gpkit.constraints.loose import Loose

Loose.reltol = 1e-4  # set the global tolerance of Loose
x = Variable('x')
x_min = Variable('x_{min}', 1)
m = Model(x, [Loose([x >= 2], senstol=1e-4),  # set the specific tolerance
              x >= x_min])
m.solve(verbosity=0)  # prints warning

Substitutions

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

Substituting into Posynomials, NomialArrays, and GPs

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

"Example substitution; 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.key: 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

"Example substitution; 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 NomialArray.

Substituting with nonnumeric values

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

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

Updating ConstraintSet substitutions

ConstraintSets have a .substitutions KeyDict attribute which will be substituted before solving. This KeyDict accepts variable names, VarKeys, and Variable objects as keys, and can be updated (or deleted from) like a regular Python dictionary to change the substitutions that will be used at solve-time. If a ConstraintSet itself contains ConstraintSets, it and all its elements share pointers to the same substitutions dictionary object, so that updating any one of them will update all of them.

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 NomialArray. 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.

Freeing Fixed Variables

After creating a Model, it may be useful to “free” a fixed variable and resolve. This can be done using the command del m.substitutions["x"], where m is a Model. An example of how to do this is shown below.

"Example of freeing fixed variables"
from gpkit import Variable, Model
x = Variable("x")
y = Variable("y", 3)  # fix value to 3
m = Model(x, [x >= 1 + y, y >= 1])
# verbosity is 0 for testing's sake, no need to do that in your code!
sol = m.solve(verbosity=0)  # optimal cost is 4; y appears in sol["constants"]

assert abs(sol["cost"] - 4) <= 1e-4
assert y in sol["constants"]

del m.substitutions["y"]
sol = m.solve(verbosity=0)  # optimal cost is 2; y appears in Free Variables
assert abs(sol["cost"] - 2) <= 1e-4
assert y in sol["freevariables"]