Unknown Bug

[codykarcher]

Wish I could be more specific here, but I'm solidly outside my area of expertise:

from gpkit import VectorVariable, Variable, units, SignomialProgram, Model
import gpkit
from gpkit.constraints.set import ConstraintSet
import numpy as np

n_Kulfan = 3

A = VectorVariable([4,n_Kulfan], "A_K", "-", "Kulfan Coefficients")

c = Variable("c", 1, "-", "Chord")
v = Variable("v", 100, "-", "Velocity")
rho = Variable("\\rho", 1.23, "-", "Density")
mu = Variable("\\mu", 1e-5, "-", "Viscosity")
Re = Variable("Re", "-", "Reynolds Number")
CL = Variable("C_L", "-", "Lift Coefficient")
CD = Variable("C_D", "-", "Drag Coefficient")

def xfoil(x0):
    Ak = np.asarray(x0["A_K"])
    CD_val = np.sum(np.sum(Ak))
    return [CD >= CD_val]

class AeroModel(ConstraintSet):
    def as_gpconstr(self, x0):
        gpconstr = super(AeroModel, self).as_gpconstr(x0)
        if x0 and x0["Re"] >= 1e7:
            constraints = xfoil(x0)
        else:
            constraints = [CD >= 1. + CL**2]
        return ConstraintSet([gpconstr, constraints])

constraints = AeroModel([CL == 1.0, Re >= c*v*rho/mu, A>=.08*np.asarray([[2, 2, 2], [1, 1, 1], [1, 1, 1], [2, 2, 1]])])

objective = CD

m = Model(objective, constraints)
sp = m.sp(require_signomials=False)
sp.localsolve('cvxopt', verbosity=1)
sp.controlpanel()

Beginning signomial solve.
Solving took 4 GP solves and 0.17 seconds.
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-15-c0a2d602937d> in <module>()
     44 m = Model(objective, constraints)
     45 sp = m.sp(require_signomials=False)
---> 46 sp.localsolve('cvxopt', verbosity=1)
     47 sp.controlpanel()

/Users/codykarcher/Dropbox (MIT)/Research/gpkit/gpkit/constraints/signomial_program.pyc in localsolve(self, solver, verbosity, x0, rel_tol, iteration_limit, **kwargs)
    128         self.sens_from_gpconstr(gp.constraints,
    129                                 result["sensitivities"]["constraints"],
--> 130                                 result["sensitivities"]["constants"])
    131         self.result = result  # NOTE: SIDE EFFECTS
    132         return result

/Users/codykarcher/Dropbox (MIT)/Research/gpkit/gpkit/constraints/set.pyc in sens_from_gpconstr(self, gpapprox, gp_sens, var_senss)
    195             gpa = gpapprox[i]
    196             gp_s = gp_sens[str(gpa)]
--> 197             constr_sens[str(c)] = c.sens_from_gpconstr(gpa, gp_s, var_senss)
    198         return constr_sens
    199 

/Users/codykarcher/Dropbox (MIT)/Research/gpkit/gpkit/constraints/set.pyc in sens_from_gpconstr(self, gpapprox, gp_sens, var_senss)
    195             gpa = gpapprox[i]
    196             gp_s = gp_sens[str(gpa)]
--> 197             constr_sens[str(c)] = c.sens_from_gpconstr(gpa, gp_s, var_senss)
    198         return constr_sens
    199 

/Users/codykarcher/Dropbox (MIT)/Research/gpkit/gpkit/constraints/set.pyc in sens_from_gpconstr(self, gpapprox, gp_sens, var_senss)
    193         constr_sens = {}
    194         for i, c in enumerate(self):
--> 195             gpa = gpapprox[i]
    196             gp_s = gp_sens[str(gpa)]
    197             constr_sens[str(c)] = c.sens_from_gpconstr(gpa, gp_s, var_senss)

IndexError: list index out of range

Only thing I do know is that it has something to do with A>=.08*np.asarray([[2, 2, 2], [1, 1, 1], [1, 1, 1], [2, 2, 1]])

[bqpd]

(this doesn't fix the problem but) You should probably switch back to the master branch, and add the following method to error model instead of using require_signomials:

    def as_posyslt1(self):
        raise ValueError("AeroModel is not allowed to solve as a GP.")

[codykarcher]

So:

from gpkit import VectorVariable, Variable, units, SignomialProgram, Model
import gpkit
from gpkit.constraints.set import ConstraintSet
import numpy as np

n_Kulfan = 3

A = VectorVariable([4,n_Kulfan], "A_K", "-", "Kulfan Coefficients")

c = Variable("c", 1, "-", "Chord")
v = Variable("v", 100, "-", "Velocity")
rho = Variable("\\rho", 1.23, "-", "Density")
mu = Variable("\\mu", 1e-5, "-", "Viscosity")
Re = Variable("Re", "-", "Reynolds Number")
CL = Variable("C_L", "-", "Lift Coefficient")
CD = Variable("C_D", "-", "Drag Coefficient")

def xfoil(x0):
    Ak = np.asarray(x0["A_K"])
    CD_val = np.sum(np.sum(Ak))
    return [CD >= CD_val]

class AeroModel(ConstraintSet):
    def as_gpconstr(self, x0):
        gpconstr = super(AeroModel, self).as_gpconstr(x0)
        if x0 and x0["Re"] >= 1e7:
            constraints = xfoil(x0)
        else:
            constraints = [CD >= 1. + CL**2]
        return ConstraintSet([gpconstr, constraints])
    def as_posyslt1(self):
        raise ValueError("AeroModel is not allowed to solve as a GP.")

constraints = AeroModel([CL == 1.0, Re >= c*v*rho/mu, A>=.08*np.asarray([[2, 2, 2], [1, 1, 1], [1, 1, 1], [2, 2, 1]])])

objective = CD

m = Model(objective, constraints)
sp.localsolve('cvxopt', verbosity=1)
sp.controlpanel()

Is that right?

[bqpd]

Almost; the last two lines are unnecessary and can be replaced by m.localsolve(). And by convention I generally put as_posyslt1 above as_gpconstr and put a line between method declarations.

[codykarcher]

Cool:

from gpkit import VectorVariable, Variable, units, SignomialProgram, Model
import gpkit
from gpkit.constraints.set import ConstraintSet
import numpy as np

n_Kulfan = 3

A = VectorVariable([4,n_Kulfan], "A_K", "-", "Kulfan Coefficients")

c = Variable("c", 1, "-", "Chord")
v = Variable("v", 100, "-", "Velocity")
rho = Variable("\\rho", 1.23, "-", "Density")
mu = Variable("\\mu", 1e-5, "-", "Viscosity")
Re = Variable("Re", "-", "Reynolds Number")
CL = Variable("C_L", "-", "Lift Coefficient")
CD = Variable("C_D", "-", "Drag Coefficient")

def xfoil(x0):
    Ak = np.asarray(x0["A_K"])
    CD_val = np.sum(np.sum(Ak))
    return [CD >= CD_val]

class AeroModel(ConstraintSet):
    def as_posyslt1(self):
        raise ValueError("AeroModel is not allowed to solve as a GP.")

    def as_gpconstr(self, x0):
        gpconstr = super(AeroModel, self).as_gpconstr(x0)
        if x0 and x0["Re"] >= 1e7:
            constraints = xfoil(x0)
        else:
            constraints = [CD >= 1. + CL**2]
        return ConstraintSet([gpconstr, constraints])

constraints = AeroModel([CL == 1.0, Re >= c*v*rho/mu, A>=.08*np.asarray([[2, 2, 2], [1, 1, 1], [1, 1, 1], [2, 2, 1]])])

objective = CD

m = Model(objective, constraints)
m.localsolve('cvxopt', verbosity=1)

Any thoughts on the original bug?

[bqpd]

Yup; instead of adding constraints, replace a default constraint, as below:

from gpkit import VectorVariable, Variable, units, SignomialProgram, Model
import gpkit
from gpkit.constraints.set import ConstraintSet
import numpy as np

n_Kulfan = 3

A = VectorVariable([4,n_Kulfan], "A_K", "-", "Kulfan Coefficients")

c = Variable("c", 1, "-", "Chord")
v = Variable("v", 100, "-", "Velocity")
rho = Variable("\\rho", 1.23, "-", "Density")
mu = Variable("\\mu", 1e-5, "-", "Viscosity")
Re = Variable("Re", "-", "Reynolds Number")
CL = Variable("C_L", "-", "Lift Coefficient")
CD = Variable("C_D", "-", "Drag Coefficient")

def xfoil(x0):
    Ak = np.asarray(x0["A_K"])
    CD_val = np.sum(np.sum(Ak))
    return CD >= CD_val

class AeroModel(ConstraintSet):
    def as_posyslt1(self):
        raise ValueError("AeroModel is not allowed to solve as a GP.")

    def as_gpconstr(self, x0):
        if x0 and x0["Re"] >= 1e7:
            self[0] = xfoil(x0)
        else:
            self[0] = CD >= 1. + CL**2
        return super(AeroModel, self).as_gpconstr(x0)

constraints = AeroModel([CD >= 1. + CL**2, CL == 1.0, Re >= c*v*rho/mu, 
                         A >= .08*np.asarray([[2, 2, 2], [1, 1, 1], [1, 1, 1], [2, 2, 1]])])

objective = CD

m = Model(objective, constraints)
_ = m.localsolve()

[codykarcher]

So how does this work if Re<1e7? It seems like I've lost the second condition...

[codykarcher]

Wait, what is self[0]? I've never seen that before...

[bqpd]

Updated with the second condition, good call. self[0] is the constraint modifying its own first constraint.

[codykarcher]

This should be renamed. I'm not sure what though @bqpd @whoburg