x0 dict changes values with iteration

[codykarcher]

It appears the x0 dictionary is changing with iteration... Is it supposed to do this?

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

# ========================================================================

class XFOIL():
    def __init__(self):
        self.pathname = "/Users/codykarcher/Xfoil/bin/./xfoil"

    def cd_model(self, x0):
        polarfile = 'polars.txt'

        # sort by alpha
        data = np.genfromtxt(polarfile, skip_header=12)
        data2 = data[data[:, 0].argsort()]

        alpha   = data2[:, 0]
        cl      = data2[:, 1]
        cd      = data2[:, 2]
        cdp     = data2[:, 3]
        cm      = data2[:, 4]
        xtr_top = data2[:, 5]
        xtr_bot = data2[:, 6]

        CD_val = np.interp(x0['C_L'], cl, cd)

        return CD >= CD_val

# ========================================================================

class CDModel(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:
            xfoil=XFOIL()
            self[0] = xfoil.cd_model(x0)
        else:
            self[0] = CD >= 1. + CL**2
        return super(CDModel, self).as_gpconstr(x0)

# ========================================================================

n_Kulfan = 3
n_afl_pts = 100
max_iter = 100

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

constraints = [CL == 1.0, 
               Re == c*v*rho/mu, 
               A  >= [[.2, .2, .2], 
                      [.1, .1, .1], 
                      [.1, .1, .1], 
                      [.2, .2, .1]],
               A <=1 
              ]

constraints = CDModel(constraints)

objective = CD

m = Model(objective, constraints)
sol=m.localsolve('cvxopt', verbosity=0)
m.controlpanel()

C_L has no lower bound
C_D has no upper bound
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-22-e3cab77f0542> in <module>()
     74 
     75 m = Model(objective, constraints)
---> 76 sol=m.localsolve('cvxopt', verbosity=0)
     77 m.controlpanel()

/Users/codykarcher/Dropbox (MIT)/Research/gpkit/gpkit/constraints/prog_factories.pyc in solvefn(self, solver, verbosity, skipsweepfailures, *args, **kwargs)
     78         else:
     79             self.program, solvefn = genfunction(self, verbosity-1)
---> 80             result = solvefn(solver, verbosity-1, *args, **kwargs)
     81             solution.append(result)
     82         solution.program = self.program

/Users/codykarcher/Dropbox (MIT)/Research/gpkit/gpkit/constraints/signomial_program.pyc in localsolve(self, solver, verbosity, x0, rel_tol, iteration_limit, **kwargs)
    103     appear to be converging, you may wish to increase the iteration limit by
    104     calling .localsolve(..., iteration_limit=NEWLIMIT).""" % len(self.gps))
--> 105             gp = self.gp(x0, verbosity=verbosity-1)
    106             self.gps.append(gp)  # NOTE: SIDE EFFECTS
    107             try:

/Users/codykarcher/Dropbox (MIT)/Research/gpkit/gpkit/constraints/signomial_program.pyc in gp(self, x0, verbosity)
    134     def gp(self, x0=None, verbosity=1):
    135         """Get a GP approximation of this SP at x0"""
--> 136         return GeometricProgram(self.cost, self.as_gpconstr(x0),
    137                                 self.substitutions, verbosity=verbosity)

/Users/codykarcher/Dropbox (MIT)/Research/gpkit/gpkit/constraints/set.pyc in as_gpconstr(self, x0)
    167 
    168         When x0 is none, may return a default guess."""
--> 169         cs = ConstraintSet([constr.as_gpconstr(x0) for constr in self])
    170         cs.substitutions.update(self.substitutions)
    171         return cs

/Users/codykarcher/Dropbox (MIT)/Research/gpkit/gpkit/constraints/set.pyc in as_gpconstr(self, x0)
    167 
    168         When x0 is none, may return a default guess."""
--> 169         cs = ConstraintSet([constr.as_gpconstr(x0) for constr in self])
    170         cs.substitutions.update(self.substitutions)
    171         return cs

<ipython-input-22-e3cab77f0542> in as_gpconstr(self, x0)
     39         if x0 and x0["Re"] >= 1e7:
     40             xfoil=XFOIL()
---> 41             self[0] = xfoil.cd_model(x0)
     42         else:
     43             self[0] = CD >= 1. + CL**2

<ipython-input-22-e3cab77f0542> in cd_model(self, x0)
     26         xtr_bot = data2[:, 6]
     27 
---> 28         CD_val = np.interp(x0['C_L'], cl, cd)
     29 
     30         return CD >= CD_val

/Users/codykarcher/Dropbox (MIT)/Research/gpkit/gpkit/keydict.pyc in __getitem__(self, key)
     94         keys = self.keymap[key]
     95         if not keys:
---> 96             raise KeyError("%s was not found." % key)
     97         values = []
     98         for key in keys:

KeyError: 'C_L was not found.'

[bqpd]

You're overwriting the first constraint with another, such that CL is no longer in the problem.

[bqpd]

Change the constraints to add your base constraint:

constraints = [CD >= 1. + CL**2,
               CL == 1.0, 
               Re == c*v*rho/mu, 

[codykarcher]

But I have CL==1, so isn't that keeping it in the problem? Or is something overwriting all my constraints?

[codykarcher]

Oh, is this the self[0] thing?

[bqpd]

Yup! The constraintset is changing its first element

[codykarcher]

Ah, it would be very ideal if it didn't have to overwrite constraints, but rather append. Especially once I get multiple high fidelity models in play... Testing now.

[bqpd]

Mm, it can append by overwriting a constraint with a ConstraintSet. From the solver's point of view, it's nice if the number of constraints contained is the same both during solving and during sensitivity generation; that's what was causing the problem earlier.

[codykarcher]

Ok... We maybe should talk about this. Actually, are you free at all tomorrow? I'm actually starting to optimize airfoils now, but I'm getting some odd behavior...

[bqpd]

Sure, taking this offline.