convexengineering / gpkit

Geometric programming for engineers
http://gpkit.readthedocs.org
MIT License
206 stars 40 forks source link

Returning Constraint in List Throwing AttributeError #615

Closed codykarcher closed 8 years ago

codykarcher commented 8 years ago

Sorry, I tried to get this down to a MWE, but can't reproduce the error:

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 Kulfan():
    def __init__(self, upper=[], lower=[]):
        self.upper = upper
        self.lower = lower

    def generate_datfile(self, afl_file='afl.dat', TE_shift = 0.0, TE_gap = 0.0005, N1=0.5, N2 = 1.0):
        ls_res = subprocess.check_output(["ls -a"], shell = True)
        ls = ls_res.split()

        if afl_file in ls_res:
            subprocess.call("rm " + afl_file, shell = True)

        Kulfan_upper = self.upper
        Kulfan_lower = self.lower
        n = len(Kulfan_upper) - 1
        n_pts = n_afl_pts
        ang = np.linspace(0,np.pi,n_pts)
        psi = (np.cos(ang)-1)/-2.
        zeta_TE = TE_shift
        zetau_T = TE_gap/2.
        zetal_T = -TE_gap/2.
        C = psi**(N1)*(1-psi)**(N2)
        S = np.zeros([n+1, psi.size])
        Su_temp = np.zeros([n+1, psi.size])
        zetau_temp = np.zeros([n+1, psi.size])
        Sl_temp = np.zeros([n+1, psi.size])
        zetal_temp = np.zeros([n+1, psi.size])

        for i in range(0,n+1):
            S[i,:]=math.factorial(n)/(math.factorial(i)*math.factorial(n-i))  *  psi**i  *  (1-psi)**(n-i)
            Su_temp[i,:] = Kulfan_upper[i]*S[i,:]
            zetau_temp[i,:] = C*Su_temp[i,:]+psi*zetau_T

            Sl_temp[i,:] = Kulfan_lower[i]*S[i,:]
            zetal_temp[i,:] = C*Sl_temp[i,:]+psi*zetal_T

        Su = np.sum(Su_temp, axis=0)
        Sl = np.sum(Sl_temp, axis=0)

        zetau = psi**(N1)*(1-psi)**(N2)*Su+psi*zetau_T
        zetal = psi**(N1)*(1-psi)**(N2)*Sl+psi*zetal_T

        camber_line = (zetau+zetal)/2.

        col1 = list(reversed(psi.tolist()))
        col2 = list(reversed(zetau.tolist()))
        col1.extend(psi[1:].tolist())
        col2.extend(zetal[1:].tolist())
        x_ud = col1
        y_ud = col2
        arr_file = np.asarray([col1, col2])
        np.savetxt(afl_file,arr_file.T)

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

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

    def cd_model(self, x0):
        Ak = x0["A_K"]
        afl_file='afl.dat'
        polarfile = 'polars.txt'

        kulfan=Kulfan(Ak[0,:]-Ak[1,:], Ak[2,:]-Ak[3,:])
        kulfan.generate_datfile(afl_file=afl_file)

        proc = subprocess.Popen([self.pathname], 
                                         stdout=subprocess.PIPE, stdin=subprocess.PIPE)
    #     proc.stdin.write('naca2412 \n' + 
        proc.stdin.write('load ' + afl_file + ' \n' + 
                         'afl \n' + 
                         'oper \n' +
                         "iter %d\n" %(max_iter)+
                         'visc \n' +
                         "%.0f \n" %(x0["Re"]) +
                         'pacc \n' +
                         '\n' +
                         '\n' +
                         'aseq 0 12 1 \n' +
                         'pwrt \n' +
                         polarfile +
                         '\n' +
                         'q')
        stdout_val = proc.communicate()[0]
#         print stdout_val
        proc.stdin.close()

        # sort by alpha
        try:
            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)
        except:
            CD_val = 1.0

        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
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-10-6ef952227131> in <module>()
    158 
    159 m = Model(objective, constraints)
--> 160 sol=m.localsolve('cvxopt', verbosity=0)
    161 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-10-6ef952227131> in as_gpconstr(self, x0)
    126         else:
    127             self[0] = CD >= 1. + CL**2
--> 128         return super(CDModel, self).as_gpconstr(x0)
    129 
    130 # ========================================================================

/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

AttributeError: 'list' object has no attribute 'as_gpconstr'
codykarcher commented 8 years ago

This works fine:

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

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

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

    def cd_model(self, x0):
        Ak = np.asarray(x0["A_K"])
        CD_val = np.sum(np.sum(Ak))
        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=1)
m.controlpanel()
codykarcher commented 8 years ago

I'm guessing it doesn't like something about this:

        # sort by alpha
        try:
            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)
        except:
            CD_val = 1.0

        return [CD >= CD_val]
codykarcher commented 8 years ago

Yea, this has the same error as the original post:

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
        try:
            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)
        except:
            CD_val = 1.0

        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()
bqpd commented 8 years ago

It's overwriting the first constraint every time it runs through as_gpcontsr

bqpd commented 8 years ago

And you should return the constraint return CD >= CD_val, not the constraints in a list (if you're only overwriting self[0]). That's probably the cause of the error?

bqpd commented 8 years ago

edited but not tested:

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
        try:
            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)
        except:
            CD_val = 1.0

        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 = [CD >= 1. + CL**2,
               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()
codykarcher commented 8 years ago

Worked thanks. Sorry again for the code dump...

codykarcher commented 8 years ago

Wait, what if I want to pass out two constraints from the model, so like:

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

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

        # sort by alpha
        try:
            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]

        return CD >= CD_val, CL>=CL_val

Or do I need a list now?

bqpd commented 8 years ago

You should return a ConstraintSet: ConstraintSet([CD >= CD_val, CL>=CL_val])

codykarcher commented 8 years ago

Ah, cool