Notebook

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