In [2]:

    

import gpkit
import gpkit.interactive
gpkit.interactive.init_printing()

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt


    

In [3]:

    

var = gpkit.Variable
vec = gpkit.VectorVariable

dCpdy = var("\\frac{dC_P}{dy}")
a = var("a")
b = var("b")
lam = var("\\lambda")
s = var("s")
xi = var("\\xi")
F = var("F")
y = var("y")
eps = var("\\epsilon")


    

In [4]:

    

gp = gpkit.GP(1/dCpdy,
              [1 >= a + b,
               2*a*b >= s*lam*y + s*xi,
               xi**2 >= (lam*y)**2 + 4*a*b,
               s*b >= dCpdy/(8*lam*y**2*F) + eps*a*b],
              {F: 1, y: 0.75, eps: 1.0/30.0})
gp


    

    
Out[4]:





$$\begin{array}[ll]
\text{}
\text{minimize}
    & \frac{1}{\frac{dC_P}{dy}} \\
\text{subject to}
    & 1 \geq a + b \\
    & 2a b \geq \lambda s y + s \xi \\
    & \xi^{2} \geq 4a b + \lambda^{2} y^{2} \\
    & b s \geq 0.125\frac{\frac{dC_P}{dy}}{\lambda F y^{2}} + a b \epsilon \\
\text{substituting}
    & F = 1 \\
    & \epsilon = 0.03333 \\
    & y = 0.75 \\
\end{array}$$

In [5]:

    

print gp.solve().table()


    

    




Using solver 'cvxopt'
Solving took 0.0341 seconds
    1.3203 : Cost
           | Free variables
\frac{dC_P}{dy} : 0.757    [-] 
   \lambda : 2.57     [-] 
       \xi : 2.14     [-] 
         a : 0.315    [-] 
         b : 0.685    [-] 
         s : 0.106    [-] 
           |
           | Constants
         F : 1        [-] 
  \epsilon : 0.0333   [-] 
         y : 0.75     [-] 
           |
           | Constant sensitivities
         F : -1       [-] 
  \epsilon : 0.11     [-] 
         y : -1       [-] 
           |

In [6]:

    

B = 3
N = 10

f = var("f")
pi = var("\\pi", np.pi)
DCp = var("(\Delta C_P)")
xi = vec(N, "\\xi")
zeta = vec(N, "\\zeta")
y = vec(N, "y")
dy = vec(N, "(\Delta y)")
lam = var("\\lambda", ('sweep', np.linspace(0.01, 8, 20)))

constraints = [1 >= a + b,
               2*a*b >= s*lam*y + s*xi,
               xi**2 >= (lam*y)**2 + 4*a*b,
               s*b >= DCp/(8*lam*dy**2*F) + eps*a*b]

constraints += [ y[-1] + dy[-1]/2 <= 1,
                [y[i] + dy[i]/2 + dy[i+1]/2 <= y[i+1] for i in range(N-1)],
                 #B*zeta[0]/(2*pi) + dy[0]/2 <= y[0],#switch0
                 dy[0]/2 <= y[0],#switch0
                 1/y >= 1 + 2*f*s/(B*a),
                 1 >= F**1.56 + 1.655*F**5.51*f**-2.76   ]

sol = gpkit.GP(1/(N*DCp), constraints, {eps: ('sweep', [0.2, 0.1, 0.05, 0.033])}).solve('mosek')
print sol.table()


    

    




Using solver 'mosek'
Sweeping 2 variables over 80 passes
Solving took 3.87 seconds
    92.607 : Cost (average of 80 values)
           | Free variables (average)
(\Delta C_P) : 0.0047   [-] 
(\Delta y) : 0.0726   [-] 
         F : 0.859    [-] 
       \xi : 2.4      [-] 
         a : 0.298    [-] 
         b : 0.702    [-] 
         f : 2.95     [-] 
         s : 0.123    [-] 
         y : 0.363    [-] 
           |
           | Constants (average)
  \epsilon : 0.0958   [-] 
   \lambda : 4        [-] 
           |
           | Constant sensitivities (average)
  \epsilon : 0.526    [-] 
   \lambda : 0.179    [-] 
           |

In [59]:

    

B = 3
N = 20
Nlams = 20

f = var("f")
pi = var("\\pi", np.pi)
DCp = var("(\Delta C_P)")
xi = vec(N, "\\xi")
zeta = vec(N, "\\zeta")
y = vec(N, "y")
dy = vec(N, "(\Delta y)")
lam = var("\\lambda", ('sweep', np.linspace(0.01, 8, Nlams)))

constraints = [1 >= a + b,
               2*a*b >= s*lam*y + s*xi,
               xi**2 >= (lam*y)**2 + 4*a*b,
               s*b >= DCp/(8*lam*dy**2*F*N/1.7) + eps*a*b]

constraints += [ y[-1] + dy[-1]/2 <= 1,
                [y[i] + dy[i]/2 + dy[i+1]/2 <= y[i+1] for i in range(N-1)],
                 #B*zeta[0]/(2*pi) + dy[0]/2 <= y[0],#switch0
                 dy[0]/2 <= y[0],#switch0
                 1/y >= 1 + 2*f*s/(B*a),
                 1 >= F**1.56 + 1.655*F**5.51*f**-2.76   ]

sol = gpkit.GP(1/(N*DCp), constraints, {eps: ('sweep', [0.2, 0.1, 0.05, 0.033])}).solve('mosek')

plt.plot(sol(lam).reshape(4, Nlams).T, sol(N*DCp).reshape(4, Nlams).T, '--', linewidth=1.5)
plt.gca().set_color_cycle(None)
Fsub = sol(F)
lamsub = sol(lam)
epssub = sol(eps)

##

lams = []
NDCps = []
for j in range(Fsub.size):
    lam = var("\\lambda", lamsub[j])

    constraints = [1 >= a + b,
                   2*a*b >= s*lam*y + s*xi,
                   xi**2 >= (lam*y)**2 + 4*a*b,
                   s*b >= DCp/(8*lam*dy**2*F*N/1.7) + eps*a*b]

    constraints += [ y[-1] + dy[-1]/2 <= 1,
                    [y[i] + dy[i]/2 + dy[i+1]/2 <= y[i+1] for i in range(N-1)],
                     #B*zeta[0]/(2*pi) + dy[0]/2 <= y[0],#switch0
                     dy[0]/2 <= y[0],#switch0
                    ]

    sol = gpkit.GP(1/(N*DCp), constraints, {F: Fsub[j], eps: epssub[j]}).solve('mosek', printing=False)
    
    lams.append(sol(lam))
    NDCps.append(sol(N*DCp))
    
lams = np.array(lams)
NDCps = np.array(NDCps)
plt.plot(lams.reshape(4, Nlams).T, NDCps.reshape(4, Nlams).T, linewidth=1.5)
plt.plot([0, 8], [16.0/27.0, 16.0/27.0], 'k')
plt.ylim([0, 0.7])
plt.xlabel('Tip Speed Ratio (lambda)')
plt.ylabel('C_P')


    

    




Using solver 'mosek'
Sweeping 2 variables over 80 passes
Solving took 4.39 seconds






    
Out[59]:





<matplotlib.text.Text at 0x7f0be748df10>






    

In [60]:

    

from IPython import utils
from IPython.core.display import HTML
import os
def css_styling():
    """Load default custom.css file from ipython profile"""
    base = utils.path.get_ipython_dir()
    styles = "<style>\n%s\n</style>" % (open(os.path.join(base,'profile_default/static/custom/custom.css'),'r').read())
    return HTML(styles)
css_styling()


    

    
Out[60]:

In [ ]: