Design Exploration and Optimization

BEMT Helicopter Hover Example

Use "Cell > Run All" to run this notebook.

import gpkit and enable $\LaTeX$ printing



In [1]:

    
%pylab inline
%config InlineBackend.figure_format = 'retina'

import gpkit
import gpkit.interactive
gpkit.interactive.init_printing()
mon = gpkit.mon
vec = gpkit.vecmon









    



Populating the interactive namespace from numpy and matplotlib

make BEM model factory for varying weight and number of bins



In [2]:

    
# Constants #
rho = mon("\\rho", 1.23, "kg/m^3", "Density of Air")
W   = mon("W", "N", "Weight of Vehicle")
xi =  vec(2,"\\xi", "N", "Thrust per bin")

# Scalars #
A     = mon("A", "m^2", "Disk Area")
Omega = mon("\\Omega", "rpm", "Rotor RPM")
R     = mon("R", "m", "Rotor Radius")
P     = mon("P", "W", "Induced Power")

# Vectors #
r   = vec(2, "r", "-", "Non-dimensional Radius")
dr  = vec(2, "\\Delta r", "-", "Non-dimensional Radius Step")
Vi  = vec(2, "V_i", "m/s", "Induced Velocities")
dCT = vec(2, "dC_T", "-", "Incremental Thrust Coefficient of Each Bin")
dCP = vec(2, "dC_P", "-", "Incremental Power Coefficient of Each Bin")
dP  = vec(2, "dP", "W", "Incremental Power of Each Bin")

def BEMT_GP(N=5, Weight=1e4):
    Weight = float(Weight)
    global xi, r, dr, Vi, dCT, dCP, dP
    
    xi =  vec(N, "\\xi", "N", "Thrust per bin")
    
    constants = {W: Weight,
                 xi: Weight/float(N) * ones(N),
                 #xi: 2*Weight/float(N)**2 * array(range(1,N+1))
                 }
    
    # Vector reasignment #
    r   = vec(N, "r", "-", "Non-dimensional Radius")
    dr  = vec(N, "\\Delta r", "-", "Non-dimensional Radius Step")
    Vi  = vec(N, "V_i", "m/s", "Induced Velocities")
    dCT = vec(N, "dC_T", "-", "Incremental Thrust Coefficient of Each Bin")
    dCP = vec(N, "dC_P", "-", "Incremental Power Coefficient of Each Bin")
    dP  = vec(N, "dP", "W", "Incremental Power of Each Bin")

    phys_constraints = (A == pi*R**2,
                        R <= 8 * gpkit.units.m,
                        Omega <= 280 * gpkit.units.rpm )

    r_constraints=([r[j] >= dr[:j].sum() for j in range(1,N)],
                   1 >= dr.sum(),
                   r[0] == dr[0]/2,
                   [r[j] >= r[j-1] + .5*dr[j-1] + .5*dr[j] for j in range(1,N)],
                   r[-1] <= 1)

    main_constraints = (xi == rho*A*(Omega*R)**2*dCT,
                        0.25 == Vi**2*r*dr/(dCT*(Omega*R)**2),
                        0.25 == Vi**3*r*dr/(dCP*(Omega*R)**3),                  
                        dP == rho*A*(Omega*R)**3*dCP,
                        P >= dP.sum())
    objective = P
    eqns = phys_constraints + r_constraints + main_constraints
    return gpkit.GP(objective, eqns, constants)

gp = BEMT_GP()

show gp



In [3]:

    
gp









    Out[3]:





$$\begin{array}[ll]
\text{}
\text{minimize}
    & P\mathrm{\left[ W \right]} \\
\text{subject to}
    & A = 3.142R^{2} \\
    & R \leq 8\mathrm{\left[ m \right]} \\
    & \Omega \leq 280\mathrm{\left[ rpm \right]} \\
    & {r}_{1} \geq {\Delta r}_{0} \\
    & {r}_{2} \geq {\Delta r}_{0} + {\Delta r}_{1} \\
    & {r}_{3} \geq {\Delta r}_{0} + {\Delta r}_{1} + {\Delta r}_{2} \\
    & {r}_{4} \geq {\Delta r}_{0} + {\Delta r}_{1} + {\Delta r}_{2} + {\Delta r}_{3} \\
    & 1 \geq {\Delta r}_{0} + {\Delta r}_{1} + {\Delta r}_{2} + {\Delta r}_{3} + {\Delta r}_{4} \\
    & {r}_{0} = 0.5{\Delta r}_{0} \\
    & {r}_{1} \geq 0.5{\Delta r}_{0} + 0.5{\Delta r}_{1} + {r}_{0} \\
    & {r}_{2} \geq 0.5{\Delta r}_{1} + 0.5{\Delta r}_{2} + {r}_{1} \\
    & {r}_{3} \geq 0.5{\Delta r}_{2} + 0.5{\Delta r}_{3} + {r}_{2} \\
    & {r}_{4} \geq 0.5{\Delta r}_{3} + 0.5{\Delta r}_{4} + {r}_{3} \\
    & 1 \geq {r}_{4} \\
    & {\xi}_{0} = {dC_T}_{0} \Omega^{2} A R^{2} \rho \\
    & {\xi}_{1} = {dC_T}_{1} \Omega^{2} A R^{2} \rho \\
    & {\xi}_{2} = {dC_T}_{2} \Omega^{2} A R^{2} \rho \\
    & {\xi}_{3} = {dC_T}_{3} \Omega^{2} A R^{2} \rho \\
    & {\xi}_{4} = {dC_T}_{4} \Omega^{2} A R^{2} \rho \\
    & 0.25 = \frac{{\Delta r}_{0} {V_i}_{0}^{2} {r}_{0}}{{dC_T}_{0} \Omega^{2} R^{2}} \\
    & 0.25 = \frac{{\Delta r}_{1} {V_i}_{1}^{2} {r}_{1}}{{dC_T}_{1} \Omega^{2} R^{2}} \\
    & 0.25 = \frac{{r}_{2} {\Delta r}_{2} {V_i}_{2}^{2}}{{dC_T}_{2} \Omega^{2} R^{2}} \\
    & 0.25 = \frac{{r}_{3} {\Delta r}_{3} {V_i}_{3}^{2}}{{dC_T}_{3} \Omega^{2} R^{2}} \\
    & 0.25 = \frac{{r}_{4} {\Delta r}_{4} {V_i}_{4}^{2}}{{dC_T}_{4} \Omega^{2} R^{2}} \\
    & 0.25 = \frac{{\Delta r}_{0} {V_i}_{0}^{3} {r}_{0}}{{dC_P}_{0} \Omega^{3} R^{3}} \\
    & 0.25 = \frac{{\Delta r}_{1} {V_i}_{1}^{3} {r}_{1}}{{dC_P}_{1} \Omega^{3} R^{3}} \\
    & 0.25 = \frac{{r}_{2} {\Delta r}_{2} {V_i}_{2}^{3}}{{dC_P}_{2} \Omega^{3} R^{3}} \\
    & 0.25 = \frac{{r}_{3} {\Delta r}_{3} {V_i}_{3}^{3}}{{dC_P}_{3} \Omega^{3} R^{3}} \\
    & 0.25 = \frac{{r}_{4} {\Delta r}_{4} {V_i}_{4}^{3}}{{dC_P}_{4} \Omega^{3} R^{3}} \\
    & {dP}_{0} = {dC_P}_{0} \Omega^{3} A R^{3} \rho \\
    & {dP}_{1} = {dC_P}_{1} \Omega^{3} A R^{3} \rho \\
    & {dP}_{2} = {dC_P}_{2} \Omega^{3} A R^{3} \rho \\
    & {dP}_{3} = {dC_P}_{3} \Omega^{3} A R^{3} \rho \\
    & {dP}_{4} = {dC_P}_{4} \Omega^{3} A R^{3} \rho \\
    & P \geq {dP}_{0} + {dP}_{1} + {dP}_{2} + {dP}_{3} + {dP}_{4} \\
\text{substituting}
    & \rho = 1.23 \\
    & {\xi}_{0} = 2000 \\
    & {\xi}_{1} = 2000 \\
    & {\xi}_{2} = 2000 \\
    & {\xi}_{3} = 2000 \\
    & {\xi}_{4} = 2000 \\
\end{array}$$

solve for local monomial model



In [4]:

    
gp.solve()["local_model"][0]









    



Using solver 'cvxopt'
Solving took 0.235 seconds






    Out[4]:





$$0.5404\frac{{\xi}_{0}^{0.32} {\xi}_{1}^{0.3} {\xi}_{2}^{0.3} {\xi}_{3}^{0.3} {\xi}_{4}^{0.27}}{\rho^{0.5}}$$



In [5]:

    
from IPython.display import Math

Math('0.5404\\frac{\\left({\\xi}_{1}^{0.0144} {\\xi}_{5}^{-0.0318} \\prod_{1}^5 {\\xi}_{i}^{0.3035} \\right)}{\\rho^{0.5}}')









    Out[5]:





$$0.5404\frac{\left({\xi}_{1}^{0.0144} {\xi}_{5}^{-0.0318} \prod_{1}^5 {\xi}_{i}^{0.3035} \right)}{\rho^{0.5}}$$

vectorized solution



In [6]:

    
gp.solution["variables"]









    Out[6]:





{V_i: array([[ 4.61764043,  4.40875243,  4.40875201,  4.4087519 ,  3.94733624]]),
 \Delta r: array([[ 0.43547418,  0.17668997,  0.13930325,  0.11877681,  0.12975581]]),
 dP: array([[ 9235.28083441,  8817.50482631,  8817.50398871,  8817.50377994,
          7894.67246202]]),
 r: array([[ 0.21773709,  0.58869686,  0.74669362,  0.87573369,  1.        ]]),
 dC_T: array([[ 0.00373315,  0.00373315,  0.00373315,  0.00373315,  0.00373315]]),
 dC_P: array([[ 0.00037037,  0.00035362,  0.00035362,  0.00035362,  0.00031661]]),
 \xi: array([[ 2000,  2000,  2000,  2000,  2000]]),
 P: array([ 43582.46513974]),
 \Omega: array([ 55.55730973]),
 A: array([ 201.06192911]),
 R: array([ 7.99999999]),
 \rho: array([ 1.23])}



In [7]:

    
gp.solution["sensitivities"]["variables"]["\\xi"]









    Out[7]:





array([[ 0.31785537,  0.30347657,  0.30347654,  0.30347653,  0.27171497]])



In [8]:

    
p = xi/rho
p









    Out[8]:





$$\begin{bmatrix}\frac{{\xi}_{0}}{\rho}\mathrm{\left[ \tfrac{N\cdotm^{3}}{kg} \right]} & \frac{{\xi}_{1}}{\rho}\mathrm{\left[ \tfrac{N\cdotm^{3}}{kg} \right]} & \frac{{\xi}_{2}}{\rho}\mathrm{\left[ \tfrac{N\cdotm^{3}}{kg} \right]} & \frac{{\xi}_{3}}{\rho}\mathrm{\left[ \tfrac{N\cdotm^{3}}{kg} \right]} & \frac{{\xi}_{4}}{\rho}\mathrm{\left[ \tfrac{N\cdotm^{3}}{kg} \right]}\end{bmatrix}$$

show effect of weight and bin number on thrust disitribution



In [10]:

    
def solve_and_plot_Vi(bins, weight):
    gp = BEMT_GP(bins, weight)
    sol = gp.solve(solver='mosek', printing=False)

    figure(figsize=(12,6))
    plot(sol(R*r), sol(xi/(R*dr)), 'r.')
    ylim(0, weight/4)
    ylabel('Thrust density [N/m]')
    xlabel('Radial position [m]')
    show()

from IPython.html.widgets import interactive

interactive(solve_and_plot_Vi, bins=(1,20), weight=(1e3, 1e5, 1e3))