convexengineering/gpkit

Geometric programming for engineers

gpkit's Issues

More mathematical 'nomials

Discussion from 93787aa

• Monomial inherits from Posynomial.
• Posynomials do not store their data as lists of Monomial terms. Instead,
• The data for a posynomial is a sparse matrix A (columns indexed by variable; rows indexed by term #). Many ways to do this, for example, we could keep the exps attribute a dictionary keyed by variable name, but make the value a list of (term#, exponent) tuples.
  • Example: the posynomial x + x*y + x*z^1.5 would be stored as {'x': [(0, 1), (1, 1), (2, 1)], 'y': [(1, 1)], 'z': [(2, 1.5)]}

fit_profile_drag.m Fails to Converge

There are a number of issues here. The first line references a non-existent director and likely should read: addpath('../convexfitting') However, even when this is corrected, the code fails to converge, getting stuck in line 7.

2D GP interpolations are seemingly not convex

Builds on #33. Could look something like:

sol = gp.splinesol({'V': np.linspace(45, 55, 50), 'V_min': np.linspace(20, 25, 50)})

Greatly reduces the number of optimizations that need to be done for a given plot.

Document install instructions

From notes taken today

Variable naming conventions

For example, in a drag model, should I create a variable D_fuselage? Or should I just put the (monomial) fuselage drag expression inside a larger posynomial? Do we offer a way for users to name terms in a posynomial? What is our recommendation for new users on best practices?

@bqpd says this is related to variable namespaces and locality.

build_solvers.py references my specific home directory

"/Users/whoburg/mosek"

Test install on Windows, OSX, and Linux

Including

• symlinking mosek and mskexpopt onto system path
• Installing dependencies (ctypesgen?)
• build instructions

etc...

Dynamically sweep Pareto frontier

Building on #30 and #34, the really cool trick would be dynamically picking local grid resolution based on the spine's smoothness after lower-resolution passes.

With a model for edge probabilities, we could target additional optimizations to minimize our uncertainty, until the estimated chance of there being a deviation greater than a defined threshold (which should have a default value based on the constraints) is below another threshold (like 0.1%).

Eyeballing the aircraft design example, I'd say that all of those graphs could be done with the right 10 points, which would only take a third of second. It might look something like this:

# skipping most declarations and imports

constants = {
    'CDA0': 0.03062702,   # [m^2] fuselage drag area
    'rho': 1.23,          # [kg/m^3] density of air
    'mu': 1.78e-5,        # [kg/(m*s)] viscosity of air
    'S_wet_ratio': 2.05,  # wetted area ratio
    'k': 1.2,             # [-] form factor
    'e': 0.96,            # [-] Oswald efficiency factor
    'W_0': 4940,          # [N] aircraft weight excluding wing
    'N_ult': 2.5,         # [-] ultimate load factor
    'V': 50,              # [m/s] cruising speed
    'V_min': 25,          # [m/s] takeoff speed
}
gpkit.monify_up(globals(), constants)

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
                ], constants=constants)

sol_i = gp.spline_solve_dynamic({'V':(45, 55, 30), 'V_min':(25, 30, 30)}) 
plt.contour(sol_i['V'], sol_i['V_min'], sol_i['C_D'])

t_GP failing (CVXopt is not properly returning the dual solution)

Very confused about why this is suddenly happening, given that unit tests were recently working, and recent commits do not include any major changes.

======================================================================
ERROR: test_trivial_gp (t_geometric_program.t_GP_cvxopt)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/whoburg/MIT/dev/gpkit/gpkit/tests/t_geometric_program.py", line 17, in test_trivial_gp
    sol = prob.solve()
  File "gpkit/geometric_program.py", line 175, in solve
    self.sensitivities = self._sensitivities(result)
  File "gpkit/geometric_program.py", line 186, in _sensitivities
    for (var, locs) in self.unsubbed.var_locs.iteritems()}
  File "gpkit/geometric_program.py", line 186, in <dictcomp>
    for (var, locs) in self.unsubbed.var_locs.iteritems()}
IndexError: index out of range

----------------------------------------------------------------------
Ran 21 tests in 0.211s

FAILED (errors=1)

Installation Instructions Links

The links to Anaconda and Mosek are broken

Units for sweep

Entering ticket from post-it note

Generalized Geometric Programs

Do we want to support min(p,p), max(p, p) and raising polynomials to positive non-integers?

From the Boyd:

Support for named vector variables / tables

From #36 discussion:

A related idea that I've played with is 'tables' (needs a better name). There are some examples in the matlab code. The idea is motivated by an aircraft design example: say I have several variables: fc_vars = gpkit.monify('W CD CL V Wfuel eta_prop'). I also have a constraint set relating these variables, already defined.

Now, I want to set up a design problem where several 'flight conditions' are considered: say, top-of-climb, cruise midway, top of descent, 3g pull-up, etc. Or, perhaps I want to simultaneously consider the max-range mission and the max-endurance mission and a max-speed mission. What I'd really like is to not have to redefine subscripted variables for each flight condition for each variable. Instead, I'd like to be able to say, there's a table, whose rows are these variables, and whose columns are these flight condition names. And I'd like to be able to apply a standard model (e.g. steady level flight) to every column in the table, so that what's going on under the hood is that when the model references V, the GP knows to constrain V_max_range and V_max_endurance and V_max_speed.

Separately, one needs the ability to constrain a single element of each of these 'vector' variables.

I'm not sure whether that made any sense - hope so.

need GP class as data model

We need a class that represents a GP (or a 'GP constraint set', if we decide we don't want it to include an objective). This ticket is a placeholder to get some ideas down on paper -- it probably needs to be split into multiple sub-tickets as we flesh these issues out.

The GP class should support:

• dict-indexing of variables, so that the column index of a variable can be looked up in O(1) time
• dict-indexing of constraints (via names??), so that the row index of a constraint can be looked up in O(1) time
• Fast re-solve for varying parameters, as enabled by indexing above.
• Combining of multiple GP constraint sets into a single GP constraint set (with coupling via shared variable names)
• Sparse matrix storage (scipy.sparse.coo_matrix) of the A matrix (cvxopt calls it F)

The data model (A, c, map) should probably match MOSEK instead of CVX, for efficiency once we get mosek working. And we presumably want ways of going back-and-forth while still supporting fast re-solve.

Cleaner, more linear history

@whoburg: The current commit history is kinda messy! I rebased my working branch; it's got half the commits and none of those confusing 'merges'.

(This is entirely cosmetic, probably.)

Friendly infeasibility

When the GP is infeasible, make sure to report:

• How infeasible is it?
• Which constraints are most problematic?
• What are some constraint changes that would make it feasible?

Mosek gp solver is currently not exposed in python

Mosek contains (to my knowledge) the fastest / best gp solver available today. Amazingly, however, despite having a Python API, MOSEK does not expose their GP solver in Python, and has no plans to do so (see e-mail transcript below).

Goal for this ticket is to implement an interface (likely using ctypes) so that MSK_expoptimize can be called from Python.

-- see test_driven_development/mosek_interface_example.py
-- see gpkit/mosek_interface.py

Possible path forward, based on digging around in the mosek Python init.py:

• compile expopt (and dgopt??) into .dylib; place in main MOSEK bin
• follow SCopt recipe in MOSEK __init__.py -- search for __scopt__
• decent ctypes reference: http://www.dalkescientific.com/writings/NBN/ctypes.html

#### BEGIN E-MAIL TRANSCRIPT

Warren Hoburg, Dec 09 03:17 (CET):
Dear Mosek Support:

I'm writing to request clarification on the API's provided for geometric
programming (also known as polynomial programming or exponential
programming) in Mosek. As far as I can tell:

• In the Mosek MATLAB toolbox, there is a function mskgpopt(c, a, map,
  ...), which allows a GP to be specified and solved directly.
• The MOSEK C API documentation describes how a Geometric Program can be
  specified as either a) a separable-convex program (exponential
  optimization), or b) a dual geometric optimization.
• The python API documentation does not appear to mention Geometric
  Programming, Exponential Optimization, or Dual Geometric programming, but
  does expose the SCopt interface.

So, my question is: if I want to solve a GP using MOSEK, and I have the c,
a, and map vectors/metrices (as defined in the MOSEK MATLAB toolbox), is
there a straightforward way to do it in Python? Would I need to go through
the motions of converting my problem to a SCopt problem? Is there a
standard or preferred way?

Thank you for any help or suggestions you can provide.

Warren Hoburg
MIT

Erling D. Andersen, Dec 09 07:26 (CET):
You cannot really solve GPs using the Python API because there is no way of feeding the dual GP.
Since it requires some work for us to make that GP from Python happen and we must support it indefinitely in the future we
are not so keen on doing that unless we fairly sure we can make some money on it.

Regards

Erling

Warren Hoburg, Dec 10 07:54 (CET):
Dear Erling,

Thank you for the prompt reply. A few more questions along these lines:

a) The documentation section 8.2.1 details the exponential optimization
formulation, but warns about numerical issues with computing large values
of e^x, and suggests using the dual geometric form instead. However, the
documentation goes on to mention options for selecting -primal or -dual for
the mskexpopt solver. So, can users reliably use the mskexpopt C api to
solve a GP (via the default dual option), or must one heed the warning and
(outside of MOSEK) convert the problem to dual form before solving it via
dgopt?

b) Wanting to solve GP's from Python, is there any reason I should not
write a simple python wrapper that calls the MOSEK C API? Are the any
difficulties you anticipate with this approach?

c) Out of curiosity, how does the MATLAB toolbox mskgpopt function work?
Does it call the mskexpopt C routines and default to solving the dual
problem? Does it do some conversion to the dual problem and then call the
dgopt C routines?

Thank you again for all your help and prompt responses.

Warren

Erling D. Andersen, Dec 10 08:10 (CET):
a) Normal the dual form is the best in our experience but sometimes the primal form wins.

b) If do that by sending c, map and A to some C code then it will be very robust. You can use the code from mskexpopt to build your C code I think. I have not thought about this but that might be your best option.

c) You can look into the MATLAB code. The general convex optimizer is available n MATLAB but is not documented.
It uses

Regards

Erling

Warren Hoburg, Dec 11 07:16 (CET):
Erling,

Thank you again for your continued help.

To clarify on your response to a), do the mskexpopt C routines handle the
decision of whether to solve primal or dual? Or is it left to the user to
recognize that the dual is usually better, convert the primal problem to
dual form, and solve via dgopt?

Warren

Erling D. Andersen, Dec 11 07:56 (CET):
There is a switch that controls whether mskexpopt solves the primal or dual formulation behind the scenes.
You input the problem the same way independdent of how MOSEK solves it.

Regards

Erling

#### END E-MAIL TRANSCRIPT

legacy matlab code -- mec Constraint is Missing

This constraint class is not in the directory.

Change name of class array

In array.py, the name array conflicts with np.array. Should probably be something that better elucidates the purpose, i.e. vectorizing calls to __geq__, __leq__, etc

Internal lookup table for PosyArray

to be filled in later...

Better variable ordering

This means something to @bqpd

Units support (to prevent silly mistakes)

Using Quantities to do this would definitely have saved me some time in the past...and units are all monomials anyway, so this shouldn't be too hard to do.

Possible bug with installation check

Possible bug: I get the error below when I use

ipython -c "%run gpkit/tests/run_tests.py"

to check my installation. I have been using the code without problem for a while, but wanted to see what the installation check threw out.

======================================================================
ERROR: test_trivial_gp (t_geometric_program.t_GP_mosek)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/philippekirschen/Documents/MIT/Research/Geometric Programming/gpkit2/gpkit/tests/t_geometric_program.py", line 17, in test_trivial_gp
    sol = prob.solve()
  File "gpkit/geometric_program.py", line 69, in solve
    result = self.__run_solver()
  File "gpkit/geometric_program.py", line 144, in __run_solver
    import _mosek.expopt
  File "gpkit/_mosek/expopt.py", line 39, in <module>
    raise ImportError("Could not load MOSEK library: "+repr(e))
ImportError: Could not load MOSEK library: ImportError('No module named expopt_h',)

Latex with subscripts greater than 9

Need braces

vector variables

Need to support vector variables. Many ways to do this -- we should start discussion.

Probably best to start with an example -- @whoburg will add to test_driven_development asap.

Variable namespaces and naming conventions

Should I be able to use V for velocity in one model, and V for voltage in another? How do we combine these models? How do we keep track of variable naming in a cross-organizational optimization model?

t_GP currently failing

Support for sweep variables in GP

Add them via GP.sweep_variables, fall back to what's in GP.constants

Whether to allow == constraint for posynomials

Posynomial equality constraints (of the form 1 == posynomial, or monomial == posynomial) are not GP compatible. Nevertheless, one could imagine offering an __eq__ operator that implies something like "relax this constraint into an inequality, solve, and confirm post-solve that equality holds".

Draft-y list of pros/cons -- to be updated according to discussion:

Pro:

• Intuitively nice for things that qualitatively can be considered definitions as opposed to constraints (like Weight = wing_weight + fuselage_weight +...

Con:

• Could confuse and ultimately frustrate users new to GP who might not realize mathematical (convexity) implications of == for posynomials. (might be more clear to receive an error earlier in the modeling process).
• Behind-the-scenes magic -- maybe functions should be provided to more explicitly do what the __eq__ operator would encode.
• More complex code

Automatic sensitivity analysis

We should do an automatic/free sensitivity analysis with the dual solution. Since we have all the constants known to the GP object, it should be easy to get useful results.

Reorganization

Idea: create new top-level examples/ directory, and move ipython notebook(s) to there?

Create aircraft design example that sweeps out Pareto frontier

Needed for apples-to-apples comparison of mosek Python interface with old results from mosek MATLAB toolbox.

CVXOPT and MOSEK both necessary?

I noticed that the installation instructions suggest installing both cvxopt and mosek. I'm not familiar with their inner workings, but I was under the impression that are substitutes for one another. I've been using the code fine (I think) with mosek only.

Just wondering if either:
a) There is a good reason to install both
b) The installation instructions should suggest them as alternatives.

Standardize our branch usage (preparation for CI)

Continuous integration: Anyone have expertise?

in flatten_constr, TypeError: 'bool' object is not iterable

Traceback (most recent call last):
  File "737_example.py", line 66, in <module>
    ], substitutions=constants, solver="mosek_cli")
  File "/Users/philippekirschen/.../gpkit/geometric_program.py", line 58, in __init__
    self.constraints = tuple(flatten_constr(constraints))
  File "/Users/philippekirschen/.../gpkit/geometric_program.py", line 21, in flatten_constr
    for elel in flatten_constr(el):
  File "/Users/philippekirschen/.../gpkit/geometric_program.py", line 17, in flatten_constr
    for el in l:
TypeError: 'bool' object is not iterable

Constants

I want to be able to define a 'model' that is a set of GP constraints. Say, for example, steady-level-flight, which might be something like:

slf = constraint_set([Monomial('L') > Monomial({'rho': 1, 'V': 2, 'CL': 1, 'S': 1}, 0.5),
                      Monomial('D') > Monomial({'rho': 1, 'V': 2, 'CD': 1, 'S': 1}, 0.5)])

Now, I'd like the ability to add this constraint set to a GP, where I have a constants dict:

constants = {'rho': 1.2}

The GP should not create a new variable 'rho'; rather, it should recognize that 'rho' is a constant, and update the c vector accordingly.

We will also need the GP model to store a matrix of exponents associated with each constant in each constraint, for use in analyzing sensitivity to each constant parameter via dual variables.

build script complains about mosek_cli

# Looking for mosek_cli
#   Trying to run mskexpopt...
#     Calling 'mskexpopt'
##
### CALL BEGINS
/bin/sh: mskexpopt: command not found
### CALL ENDS
##
# Did not find mosek_cli

Iterated relaxation for integer variables

Which should have the unit "count"

More examples would be great

Talking to friends who could probably use this in their jobs but aren't used to thinking of things this way, it seems like the mose useful introduction would be many small-but-interesting examples. A standard format for them might make them easier to parse, as well.

gp.m print_report calls undefined function

Undefined function 'print_table' for input arguments of type 'cell'.

Error in gp/print_report (line 356)
print_table(vals, vs, fc);

Error in drone3 (line 116)
p.print_report(res, 'print_sensitivities', true);

Support for vector constants in GP

Automatically expand vector-valued constants (e.g. x = [1, 2, 3]) into (x0, x1, x2 = 1, 2, 3) in Monomial.sub, and in GP.constants_update

Separate latex representation and python variable name?

 In: W_rbar = Monomial('\bar{W_r}') 
     print W_rbar
Out: Monomial(\bar{W_r})

Should this be allowed?

Test_Monomial unit test currently failing

ERROR: test_init (main.Test_Monomial)

Traceback (most recent call last):
File "test_gpkit.py", line 13, in test_init
self.assertEqual(m.a, [1, -1])
AttributeError: 'Monomial' object has no attribute 'a'

Ran 8 tests in 0.002s

FAILED (errors=1)

Build script

We need something that:

1. If MOSEK is installed, attempts to:
   • Build build/mosek_h.py by running ctypesgen
   • Build the expopt shared library into build/libexpopt.so by pulling the C files, diffing them, and compiling
2. Runs tests with the various solvers
3. Writes a config file with a default solver chosen

This might be run whenever gpkit is imported, if you don't already have a config file, or (perhaps more properly) has to be run explicitly.

Upload legacy m code

Further abstract models

There's a potential for two more abstract classes:

• Things Which Can be Substituted Into (which could also be the LaTex-printing parent)
• Things Which Have exps & cs Attributes

ConstraintSet object

There was a discussion on this here: 22db4bb

Wanted:

• immutability
• indexing by variable
• names for constraints and cost functions
  • default name: during inequalities, the string representation of the monomial side
    • e.g. (D >= 0.5*rho*C_D*S*V**2) is gpkit.Constraint(D**-1 * 0.5*rho*C_D*S*V**2, 'D')
  • failing that, the posynomial's string representation

The below implements immutability:

class ConstraintSet(tuple):
    def __add__(self, other):
        return list(self) + other

    def __radd__(self, other):
        return other + list(self)

import gpkit
rho, mu, V, C_L, C_D, S, A, W, D, Re = gpkit.monify('rho mu V C_L C_D S A W D Re')

lift = 0.5*rho*C_L*S*V**2
drag = 0.5*rho*C_D*S*V**2
Re_lim = (rho/mu)*V*(S/A)**0.5

steady_level_flight = ConstraintSet(W <=lift,
                                    D >= drag,
                                    Re <= Re_lim)

assert steady_level_flight[1] == (W <= lift)
steady_level_flight[1] = D  # fails! cannot be modified
steady_level_flight + [W > 1520]  # returns a list, whose elements can be modified...

speed comparisons: MOSEK ~10% different than MATLAB and cvxopt

It seems (anecdotally) as though MOSEK might run much faster in Matlab than it does in gpkit. Let's test this, and if it's true bring a bug to MOSEK.

Missing NACA24xx_fits

This file is missing from the directory. May stem from other issues on this list.

Do automatic solution checking and validation

Very low priority:

We could create GP.solution_check(self, sol). It would validate that a particular primal/dual solution pair solves the GP. Could be useful for sanity checks and trying out 'suspect' solution codes.

This is very very very low priority.