convexengineering / gpkit

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

SP program not working with Model program #406

Closed bs1912 closed 7 years ago

bs1912 commented 8 years ago

I'm having problems with a SP program I developed with the previous GPKit version. I wonder whether someone could help me out, as simply replacing the SP function with the Model function doesn't do the trick. Below is the original code without any modifications. I appreciate the help!

Bruno Scalzo

#!/usr/bin/python

import gpkit
import numpy as np
import math as math
from gpkit import Variable, VectorVariable, GP
gpkit.disable_units()

# Constants
pi               = Variable("\\pi", np.pi, "-", "Half of the circle constant")
g                = Variable("g", 9.81, "m/s^2", "Gravitational constant")
nu               = Variable("\\nu", 16.0e-6, "m^2/s", "Kin. viscosity at sea level")
rho              = Variable("\\rho", 1.225, "kg/m^3", "Air density at sea level")
eps0             = Variable("\\epsilon_{0}", 8.85418782e-12, "m^-3kg_1s^4A^2", "permitivity of free space")

# Fuselage variables
l_f                = Variable("l_f", 1.0668, "m", "fuselage length")
t_f                = Variable("t_f", 0.1905, "m", "fuselage thickness")

# Wing variables
b                = Variable("b", "m", "wing span")
c                = Variable("c", "m", "wing chord")
nr               = Variable("nr", "-", "number of ribs")
rs               = Variable("rs", 0.15, "m", "rib spacing")
tau              = Variable("\\tau", 0.12, "-", "wing percent thickness")
cosdi            = Variable("cosdi", 0.939693, "-", "cosine of dihedral")
rt               = Variable("rt", 0.00635, "m", "rib thickness")
st               = Variable("st", 0.04, "-", "% thickness of spar normalized by chord")
cft              = Variable("cft", 1.778e-4, "m", "carbon fiber thickness")
SA_ref           = Variable("SA_{ref}", 2.05967,"-", "Wing SA ref. area")
Ale_ref          = Variable("Ale_{ref}", 0.010668,"-", "Wing LE ref. area")
Ate_ref          = Variable("Ate_{ref}", 0.006699,"-", "Wing TE ref. area")
Ar_ref           = Variable("Ar_{ref}", 0.085379,"-", "Wing rib ref. area")
ml_sd            = Variable("ml_{sd}", 0.020343,"kg/m^2", "microlite surface density")
rho_h60          = Variable("\\rho_{h60}", 35.24, "kg/m^3", "Density of high load 60 foam")
rho_f250         = Variable("\\rho_{f250}", 25, "kg/m^3", "Density of foamular250")
rho_cf           = Variable("\\rho_{cf}", 1600, "kg/m^3", "Density of carbon fibre sheet")
l_w              = Variable("l_{w}", "m", "Wing moment arm")
C_L_max          = Variable("C_L_max", 1.5312, "-", "Maximum Lift coefficient from tS12 XFOIL analysis")

# HT variables
A_ht             = Variable("A_{ht}", 5, "-", "HT aspect ratio")
A_vt             = Variable("A_{vt}", 5, "-", "VT aspect ratio")
l_ht             = Variable("l_{ht}", 1.5, "m", "HT moment arm")
l_vt             = Variable("l_{vt}", 1.5, "m", "VT moment arm")
S_ht             = Variable("S_{ht}", "-", "HT planform area")
S_vt             = Variable("S_{vt}", "-", "VT planform area")
c_ht             = Variable("c_{ht}", "m", "HT chord")
c_vt             = Variable("c_{vt}", "m", "VT chord")
b_ht             = Variable("b_{ht}", "m", "HT span")
b_vt             = Variable("b_{vt}", "m", "VT span")
Vh               = Variable("Vh", "-", 0.6, "tail volume coeff., HT")
Vv               = Variable("Vv", "-", 0.05, "tail volume coeff., VT")
tau_ht           = Variable("\\tau_{ht}", "-", "HT thickness to chord ratio")
C_L_ht           = Variable("C_{L_{ht}}", "-", "HT Lift coefficient")
Re_ht            = Variable("Re_{ht}", "-", "HT Reynolds number")

# Weight variables
W_batt           = Variable("W_{batt}", "N", "Battery Weight")
Mb               = Variable("Mb", 0.0038, "kg", "Single cell battery mass")
n_batt           = Variable("n_batt", "-", "number of batteries")
eta_pc           = Variable("\\eta_pc", 0.85, "W/W", "power conversion efficiency")
W_pc             = Variable("W_{pc}","N", "Power Converter Weight")
W_f              = Variable("W_{f}", 1.1772,"N", "Weight of fuselage")
W_w              = Variable("W_{w}","N", "Wing weight")
W_s              = Variable("W_{s}","N", "Wing weight, spar")
W_le             = Variable("W_{le}","N", "Wing weight, LE")
W_te             = Variable("W_{te}","N", "Wing weight, TE")
W_r              = Variable("W_{r}","N", "Wing weight, ribs")
W_ml             = Variable("W_{ml}","N", "Wing weight, microlite")
W_t              = Variable("W_{t}","N", "Thruster weight")
W_ts             = Variable("W_{ts}","N", "Tail spar weight")
tspm             = Variable("tspm", 0.0316964,"kg/m", "Tail spar mass per m")
sspm2            = Variable("sspm2", 0.19375, "kg/m^2", "side strut mass per m^2")
W_ht             = Variable("W_{ht}","N", "HT weight")
W_vt             = Variable("W_{vt}","N", "VT weight")
W_ss             = Variable("W_{ss}","N", "side strut weight")
A_ss             = Variable("A_{ss}", 0.15484, "m^2", "side strut area")
W                = Variable("W", "N", "Aircraft's  weight")

# Thruster variables
T                = Variable("T", "N", "Thrust force")
ndi_frac         = Variable("ndi_{frac}", 0.6, "-", "percent not-dihedral")
di_frac          = Variable("di_{frac}", 0.4, "-", "percent dihedral")
di_corr          = Variable("di_{corr}", 0.975877, "-", "percent dihedral")
P                = Variable("P", "kW", "Power required")
N                = Variable("N", "-", "Number of thrusters/stage")
K                = Variable("K", "-", "Number of stages")
c_t              = Variable("c_t", 2*0.0254,"m", "chord length of thruster")
e_d              = Variable("e_d", 0.202e-3,"m", "emitter diamaeter")
C_De             = Variable("C_De", 2,"-", "emitter drag coefficient")
A_t              = Variable("A_t", "m^2", "Thruster characteristic area")
L_t              = Variable("L_t", "m", "Thruster characteristic length")
Ke_fac           = Variable("Ke_fac", "-", "Multiplier for stage drag effects, emitter electrode")
Kg_fac           = Variable("Kg_fac", "-", "Multiplier for stage drag effects, ground electrode")
E_t              = Variable("E_t", "V/m", "Thruster electric field")
V_iw             = Variable("V_iw", "m/s", "ionic wind velocity")
G_t              = Variable("G_t", 0.05, "m", "Thruster gap distance")
T2P              = Variable("T2P", "N/kW", "Thruster thrust-to-power ratio")
V_op             = Variable("V_{op}", "kV", "Thruster operating voltage")
DeltaL           = Variable("DeltaL", "-", "non-dimensional pair spacing")
TA               = Variable("TA", "N/m^2", "thrust per unit area")
Trat             = Variable("Trat", "-", "thrust ratio")
T1               = Variable("T1", "N", "thrust of single pair")
Gain             = Variable("Gain", 200, "-", "Voltage Gain of transformer")
eta_prop         = Variable("eta_prop", "-", "Propulsive efficiency")

# Drag variables
C_Dp             = Variable("C_{D_p}", "-", "drag model parameter")
C_DA             = Variable("C_DA", "m^2", "Total vehicle drag times area")
C_D              = Variable("C_D", "-", "Total vehicle drag coefficient")
C_Dt             = Variable("C_Dt", "-", "Thruster drag coefficient")
C_Dss            = Variable("C_Dss", "-", "Side strut drag coefficient")
C_Di_w           = Variable("C_{D_iw}", "-", "induced drag coefficient, wing")
C_DA_F           = Variable("C_{D}A_{F}", "m^2", "Fuselage drag coefficient * frontal area")
C_D_ht           = Variable("C_{D_{ht}}", "-", "Horizontal Tail Drag coefficent")
C_Di_ht          = Variable("C_{D_{i_{HT}}}", "-", "induced drag coefficient, HT")
e                = Variable("e", 0.95, "-", "Wing spanwise efficiency")
Re               = Variable("Re", "-", "Reynold's number")
Re_f             = Variable("Re_{f}", "-", "Reynold's number, for fuselage (MTFLOW)")
Re_ss            = Variable("Re_{ss}", "-", "Reynold's number, for side strut")
Re_t             = Variable("Re_{t}", "-", "Reynold's number, for thruster")
D                = Variable("D", "N", "Total Drag on aircraft")
Dt               = Variable("Dt", "N", "Drag on thruster")
Dnt              = Variable("Dnt", "N", "Drag not on thruster")

# Structural variables
N_lift           = Variable("N_{lift}", "N", "Wing loading multiplier")
sig_maxh60       = Variable("\\sigma_{max,h60}", 15e3, "Pa", "Allowable stress, highload 60")
sigs_maxh60      = Variable("\\sigma_{max,shear,h60}", 400e3, "Pa", "Allowable shear stress")
sig_maxcf        = Variable("\\sigma_{max,cf}", 570e6, "Pa", "Allowable stress, carbon fibre")
sigs_maxcf       = Variable("\\sigma_{max,shear,cf}", 90e6, "Pa", "Allowable shear stress")
E_h60            = Variable("E_{h60}", 15000e3, "Pa", "Young's Modulus of high load 60 foam")
E_cf             = Variable("E_{cf}", 140e9, "Pa", "Young's Modulus of carbon fibre")
G_h60            = Variable("G_{h60}", 7.5e6, "Pa", "Shear Modulus high load 60 foam")
delta_max        = Variable("\\delta_{max}}", 0.05, "m", "Maximum tip deflection of wing")
Mom_w            = Variable("Mom_{w}", "Nm", "root bending moment of wing")
MomW_w           = Variable("MomW_{w}", "Nm", "root bending moment of wing, due to weight")
I_w              = Variable("I_{w}", "m^4", "second area moment of wing spar")
Q                = Variable("Q", "N", "shear force in wing spar")
QW               = Variable("QW", "N", "shear force in wing spar due to weight")
delw_b           = Variable("delw_{b}", "m", "deflection due to bending")
delw_s           = Variable("delw_{s}", "m", "deflection due to shear")
del_max          = Variable("del_{max}", 0.025, "-", "max % displacement allowed")
delw             = Variable("delw", "m", "total deflection of wing")

# Flight variables
V                = Variable("V", "m/s", "Flight speed")
C_L              = Variable("C_L", "-", "Wing lift coefficent")
L                = Variable("L", "N", "Total wing lift")
S                = Variable("S", "m^2", "Wing area")
A                = Variable("A", "-", "Aspect Ratio")
gam              = Variable("\\gamma","rad","Flight angle")
E                = Variable("E", "-", "Aerodynamic efficiency")

# Stability variables
xnp              = Variable("xnp","m","neutral point location")
xcg              = Variable("xcg","m","center of gravity location")
zcg              = Variable("zcg", "m", "center of gravity location in z-axis")
C_M              = Variable("C_{M}", -0.0440, "-", "Averaged Moment coefficient of tS12 airfoil")
C_M_ht           = Variable("C_{M_{HT}}", -0.0453, "-", "Averaged moment coefficient of HT")
htf              = Variable("htf", "-", "Vh*S/(l_ht*c)")

# Range equations variables
E_spec           = Variable("E_spec", "W*h/kg", "Mass specific energy content of battery")
Range            = Variable("Range", "m", "Range of UAV")
Endurance        = Variable("Endurance", "s", "Endurance of flight")
Rmin             = Variable("R_min", "m", "Minimum turn radius")

# Sweep variables
alfa1            = Variable("alfa_1", 0.4, "-", "Sweep coefficient")
alfa2            = Variable("alfa_2", 0.6667, "-", "Sweep coefficient")
alfa3            = Variable("alfa_3", 0.6667, "-", "Sweep coefficient")

# Utopian values for objective function
Range_ut         = Variable("Range_ut", 2459.09, "m", "Utopian Range")
E_ut             = Variable("E_ut", "-", 18.60, "Utopain aerodonyamic efficiency")
W_ut             = Variable("W_ut", 10.29, "N", "Utopian weight")
eta_prop_ut      = Variable("eta_prop_ut", 0.0505 ,"-", "Utopian propulsive efficiency")

gpkit.enable_signomials()

eqns     = [

#### Flight Model ####

                L == 0.5*rho*C_L*S*V**2*di_corr, # define lift
                Dnt == 0.5*rho*V**2*C_DA, # define drag, everything but thruster
                Dt >= 0.5*rho*V**2*(C_Dt*K*A_t+C_De*L_t*e_d) +
                          0.5*rho*V_iw**2*(C_Dt*Kg_fac*A_t+C_De*Ke_fac*L_t*e_d), # define drag, thruster
                D >= Dnt + Dt, # define total drag

                # Impose condition for Stall C_L vs Re relationship
                Re >= 176.8/C_L**-17.5 + 7.429e5/C_L**-2.275 - 7.425e5/C_L**-2.129,

                Re == (1/nu)*V*c, # define wing Re
                Re_f == (1/nu)*V*0.10668, # define fuselage Re (MTFLOW)
                Re_t == (1/nu)*V*2*0.0254, # define ground electrode Re
                Re_ss == (1/nu)*V*20*0.0254, # define side strut Re
                Re_ht == (1/nu)*V*c_ht, # define HT Re

                ## Steady Level Flight ##
                L >= W,
                T >= D,

                ## Climbing ##
#               T >= W*gam + D, 
#               L >= L,

                b == (S*A)**0.5, # define wing span
                c == (S/A)**0.5, # define wing chord
                C_D == D/(0.5*rho*V**2*S), # calculate vehicle drag coefficient  

                E == L/D, # Aerodynamic efficiency

#### Thrust ####

                T == (A*S)**0.5*ndi_frac*N*K/(T1*Trat), # get thrust given thruster geometry
                T1 >= 1221695966.01*V_op**-7.3747673794*1.e3/5.,
                T1 >= 4151.46592088*V_op**-3.12069777202*1.e3/5.,
                T1 >= 67112.3207544*V_op**-3.98906732575*1.e3/5.,
                Trat**54.74862 >= 0.06725*DeltaL**-13.99422 + 0.91309*DeltaL**-1.92270,
                P == T/T2P, # estimate power from T/P
                P <= 0.750, # limited by expected output from PC prototype
                E_t == V_op*1.e3/G_t, # calculate thruster electric field (1D)
                V_iw == (E_t**2*eps0/rho)**0.5, # calculate ionic wind velocity
                K >= 1,
                N >= 1,
                K <= 3,
                T2P == G_t/(2e-4*V_op),
                DeltaL*N*G_t <= t_f,
                V_op >= 15,
                V_op <= 40,
                n_batt == V_op*1000/(Gain*3.3),
                DeltaL <= 2.0,
                DeltaL >= 0.4,

                eta_prop == T*V/(P*1000),

#### Drag Model ####

                C_DA_F**2 == 4.897e-2*Re_f**-0.9811, # fuselage drag

                C_Di_w >= C_L**2/(pi*A*e) + 0.38*C_Dp*C_L**2, # wing induced drag (Schaufelle model)

                C_Di_ht >= C_L_ht**2/(pi*A*e) + 0.38*C_D_ht*C_L_ht**2, # HT induced drag (Schaufelle model)

                C_Dss == 1.33/Re_ss**0.5, # side strut drag

                1 >= (0.0309*C_L**4.9014/(Re**1.9359*C_Dp**5.826) +
                     166.9188*C_L**20.4236/Re**3.4473*C_Dp**8.1109) +
                     0.1722*Re**0.0887/(C_L**0.0722*C_Dp**0.0957) +
                     1.8286e7*C_L**0.1999/(Re**2.4252*C_Dp**2.3477), # wing profile drag

                A_t == ndi_frac*(A*S)**0.5*c_t*N, # planform area, ground electrodes
                L_t == ndi_frac*(A*S)**0.5*N, # length of electrodes per stage
                Kg_fac == 0.921594*K**1.75391, # factor to account for increasing speed between stages
                Ke_fac == Kg_fac, # fit no good, assume this as "worst case"

                # thruster drag
                C_Dt**10.5904 >= 5.6706e-9*Re_t**-1.5691 + 4750.3152*Re**-4.4195 + 7.4263e-6*Re**-2.9821,

                # vehicle drag coeff. times area, no thruster
                C_DA >= C_DA_F  + C_Dp*S + C_Di_w*S + 2*2*C_Dss*A_ss + C_Di_ht*S_ht             + C_D_ht*S_ht,

#### Weight Constraints ####

                W >=  W_f + W_w + W_batt + W_pc + W_t + W_ht + W_vt + W_ts + W_ss, # total weight
                P/eta_pc/n_batt <= 15.e-3, # limit cell power
                W_batt == g*n_batt*Mb, # battery weight
                nr == b/rs, # calculate number of ribs in wing
                W_ml == g*SA_ref*c*b*ml_sd, # weight of microlite
                W_le == g*Ale_ref*c**2*b*rho_f250, # weight of LE
                W_te == g*Ate_ref*c**2*b*rho_f250, # weight of TE
                W_r == g*Ar_ref*c**2*rt*nr*rho_f250, # weight of ribs
                W_s >= g*(c**2*tau*st)*b*rho_h60 + g*cft*(st*c)*b*rho_cf, # weight of spar
                W_w >= W_ml + W_te + W_le + W_r + W_s, # wing weight
                W_t == g*0.0134*(A*S)**0.5*ndi_frac*N*K, # thruster weight
                W_ts == g*l_ht*tspm, # tail spar weight
                W_ht == W_w*S_ht/S, # HT weight
                W_vt == W_w*S_vt/S, # VT weight
                W_pc == g*0.375, # power converter weight
                W_ss == 2*g*sspm2*A_ss, # side strut weight

#### Structural Model ####

                Mom_w == L/2*b/4, # calculate root bending moment, from lift (worst case, in flight)
                # calculate root bending moment, from weight (not in flight, worst case)
                MomW_w >= W_w/2*b/4 + W_ss/2*ndi_frac*b/2 + W_t/2*ndi_frac*b/4,
                I_w == 2*(st*c)*cft*(tau*c/2)**2, # spar moment of inertia (just cf)
                sig_maxcf >= (N_lift*Mom_w*c*tau/2)/I_w, # bending stress limit
                sig_maxcf >= (N_lift*MomW_w*c*tau/2)/I_w, # bending stress limit, weight
                Q == L/2, # calculate root shear stress
                QW >= W_w/2 + W_ss/2 + W_t/2, # calculate root shear stress due to weight
                sigs_maxh60 >= N_lift*Q/(st*tau*c**2), # shear stress limit (just foam)
                sigs_maxh60 >= N_lift*QW/(st*tau*c**2), # shear stress limit (just foam), weight
                delw_b == N_lift*L/2*(b/2)**3/(8*E_cf*I_w), # deflection due to bending
                delw_s == N_lift*Mom_w/(2*G_h60*(st*tau*c**2)), # deflection due to shear
                delw >= delw_b + delw_s, # total deflection
                delw/b <= del_max, # deflection limit

#### Turning requirements ####

                N_lift == 0.5*rho*V**2*C_L_max/(W/S), # C_L_max for tS12 = 1.5312

                Rmin == 1.45*V**2/(g*N_lift**1.299), # minimum turn radius (posynomial approximation)

                Rmin <= 50/pi, # minimum turning radius must be smaller than track field radius

#### Stability requirements ####

# get required tail lift coefficient

                # Moment arms for now
                zcg == 0.01, # must be defined (cannot = 0)
                l_w == 0.1,

                # Lift coefficient requirement
                C_L >= C_L_ht,

                # xcp_w < xcg
                C_L_ht*htf*l_ht >= l_w*C_L + c*C_M + c_ht*C_M_ht*htf + zcg*T/L,

# get drag coefficient of HT

                1 >= 0.2201*C_D_ht**-0.2762*C_L_ht**-0.078*Re_ht**-0.3878*tau_ht**0.8912
                        + 0.0046*C_D_ht**-0.2109*C_L_ht**0.472*Re_ht**0.2637*tau_ht**-0.1398
                        + 2.1243*C_D_ht**-0.4001*C_L_ht**-0.0241*Re_ht**-0.2605*tau_ht**0.1337,

                tau_ht >= 0.08,
                tau_ht <= 0.14,

#### Tail properties ####

# get HT span, chord, and area
                S_ht == Vh*S*c/l_ht,
                b_ht == (A_ht*S_ht)**0.5,
                c_ht == S_ht/b_ht,

                htf == Vh*c/(l_ht),

# get VT span, chord, and area
                S_vt == Vv*S*b/l_vt,
                b_vt == (A_vt*S_vt)**0.5,
                c_vt == S_vt/b_vt,
#### COG Calculations ####

# get neutral point

                (c/xnp) == 2.18270636553*A**-0.106088939277,

#### Range Equations ####

                E_spec == (0.33)/Mb, # Wh/kg
                Range == V*(W_batt/g)*E_spec/(P*1000)*60**2,
                Endurance == Range/V,

                ]

#### Objective function ####

Obj = alfa1*alfa2*alfa3*(W/W_ut) + alfa2*alfa3*(E_ut/E) - alfa1*alfa2*alfa3*(E_ut/E) + alfa3*(Range_ut/(Range)) - alfa2*alfa3*(Range_ut/(Range)) + (eta_prop_ut / (T*V/(P*1000))) - alfa3*(eta_prop_ut / (T*V/(P*1000)))

##### Single objective #####
sp = gpkit.SP(Obj, eqns)
sol = sp.localsolve()
#####

print sol.table()
bqpd commented 8 years ago

The problem appears to have come in with the new Signomial constraints, because 945cfd9e85d57f29934f229184a66dccbf020351 solves it (getting a faster solution, slightly better but within reltol).

bqpd commented 8 years ago

This works with the SP initial guess of "treat signomials as posynomial constraints, ignoring negative terms". The issue appears to be this constraint with large constants: Re >= 176.8/C_L**-17.5 + 7.429e5/C_L**-2.275 - 7.425e5/C_L**-2.129. This was previously evaluated as 1 >= 176.8/C_L**-17.5/Re + 7.429e5/C_L**-2.275/Re - 7.425e5/C_L**-2.1/Re; now it is treated as Re + 7.425e5/C_L**-2.129 >= 176.8/C_L**-17.5 + 7.429e5/C_L**-2.275

bqpd commented 8 years ago

Oh, the constraint C_L_ht*htf*l_ht >= l_w*C_L + c*C_M + c_ht*C_M_ht*htf + zcg*T/L becomes a Signomial after substitution (neat to see that), and also has to be changed to 0 >= l_w*C_L/(C_L_ht*htf*l_ht) + c*C_M/(C_L_ht*htf*l_ht) + c_ht*C_M_ht*htf/(C_L_ht*htf*l_ht) + zcg*T/L/(C_L_ht*htf*l_ht) - 1. The 1 >= form doesn't work properly in the conversion from a Posynomial constraint to a Signomial constraint.

bqpd commented 8 years ago

@bs1912, the following code now works in the newest development version, but is clearly less readable than the original, and that should be fixed before this issue is closed.

import gpkit
import numpy as np
import math as math
from gpkit import Variable, VectorVariable, GP
gpkit.disable_units()

# Constants
pi               = Variable("\\pi", np.pi, "-", "Half of the circle constant")
g                = Variable("g", 9.81, "m/s^2", "Gravitational constant")
nu               = Variable("\\nu", 16.0e-6, "m^2/s", "Kin. viscosity at sea level")
rho              = Variable("\\rho", 1.225, "kg/m^3", "Air density at sea level")
eps0             = Variable("\\epsilon_{0}", 8.85418782e-12, "m^-3kg_1s^4A^2", "permitivity of free space")

# Fuselage variables
l_f                = Variable("l_f", 1.0668, "m", "fuselage length")
t_f                = Variable("t_f", 0.1905, "m", "fuselage thickness")

# Wing variables
b                = Variable("b", "m", "wing span")
c                = Variable("c", "m", "wing chord")
nr               = Variable("nr", "-", "number of ribs")
rs               = Variable("rs", 0.15, "m", "rib spacing")
tau              = Variable("\\tau", 0.12, "-", "wing percent thickness")
cosdi            = Variable("cosdi", 0.939693, "-", "cosine of dihedral")
rt               = Variable("rt", 0.00635, "m", "rib thickness")
st               = Variable("st", 0.04, "-", "% thickness of spar normalized by chord")
cft              = Variable("cft", 1.778e-4, "m", "carbon fiber thickness")
SA_ref           = Variable("SA_{ref}", 2.05967,"-", "Wing SA ref. area")
Ale_ref          = Variable("Ale_{ref}", 0.010668,"-", "Wing LE ref. area")
Ate_ref          = Variable("Ate_{ref}", 0.006699,"-", "Wing TE ref. area")
Ar_ref           = Variable("Ar_{ref}", 0.085379,"-", "Wing rib ref. area")
ml_sd            = Variable("ml_{sd}", 0.020343,"kg/m^2", "microlite surface density")
rho_h60          = Variable("\\rho_{h60}", 35.24, "kg/m^3", "Density of high load 60 foam")
rho_f250         = Variable("\\rho_{f250}", 25, "kg/m^3", "Density of foamular250")
rho_cf           = Variable("\\rho_{cf}", 1600, "kg/m^3", "Density of carbon fibre sheet")
l_w              = Variable("l_{w}", "m", "Wing moment arm")
C_L_max          = Variable("C_L_max", 1.5312, "-", "Maximum Lift coefficient from tS12 XFOIL analysis")

# HT variables
A_ht             = Variable("A_{ht}", 5, "-", "HT aspect ratio")
A_vt             = Variable("A_{vt}", 5, "-", "VT aspect ratio")
l_ht             = Variable("l_{ht}", 1.5, "m", "HT moment arm")
l_vt             = Variable("l_{vt}", 1.5, "m", "VT moment arm")
S_ht             = Variable("S_{ht}", "-", "HT planform area")
S_vt             = Variable("S_{vt}", "-", "VT planform area")
c_ht             = Variable("c_{ht}", "m", "HT chord")
c_vt             = Variable("c_{vt}", "m", "VT chord")
b_ht             = Variable("b_{ht}", "m", "HT span")
b_vt             = Variable("b_{vt}", "m", "VT span")
Vh               = Variable("Vh", "-", 0.6, "tail volume coeff., HT")
Vv               = Variable("Vv", "-", 0.05, "tail volume coeff., VT")
tau_ht           = Variable("\\tau_{ht}", "-", "HT thickness to chord ratio")
C_L_ht           = Variable("C_{L_{ht}}", "-", "HT Lift coefficient")
Re_ht            = Variable("Re_{ht}", "-", "HT Reynolds number")

# Weight variables
W_batt           = Variable("W_{batt}", "N", "Battery Weight")
Mb               = Variable("Mb", 0.0038, "kg", "Single cell battery mass")
n_batt           = Variable("n_batt", "-", "number of batteries")
eta_pc           = Variable("\\eta_pc", 0.85, "W/W", "power conversion efficiency")
W_pc             = Variable("W_{pc}","N", "Power Converter Weight")
W_f              = Variable("W_{f}", 1.1772,"N", "Weight of fuselage")
W_w              = Variable("W_{w}","N", "Wing weight")
W_s              = Variable("W_{s}","N", "Wing weight, spar")
W_le             = Variable("W_{le}","N", "Wing weight, LE")
W_te             = Variable("W_{te}","N", "Wing weight, TE")
W_r              = Variable("W_{r}","N", "Wing weight, ribs")
W_ml             = Variable("W_{ml}","N", "Wing weight, microlite")
W_t              = Variable("W_{t}","N", "Thruster weight")
W_ts             = Variable("W_{ts}","N", "Tail spar weight")
tspm             = Variable("tspm", 0.0316964,"kg/m", "Tail spar mass per m")
sspm2            = Variable("sspm2", 0.19375, "kg/m^2", "side strut mass per m^2")
W_ht             = Variable("W_{ht}","N", "HT weight")
W_vt             = Variable("W_{vt}","N", "VT weight")
W_ss             = Variable("W_{ss}","N", "side strut weight")
A_ss             = Variable("A_{ss}", 0.15484, "m^2", "side strut area")
W                = Variable("W", "N", "Aircraft's  weight")

# Thruster variables
T                = Variable("T", "N", "Thrust force")
ndi_frac         = Variable("ndi_{frac}", 0.6, "-", "percent not-dihedral")
di_frac          = Variable("di_{frac}", 0.4, "-", "percent dihedral")
di_corr          = Variable("di_{corr}", 0.975877, "-", "percent dihedral")
P                = Variable("P", "kW", "Power required")
N                = Variable("N", "-", "Number of thrusters/stage")
K                = Variable("K", "-", "Number of stages")
c_t              = Variable("c_t", 2*0.0254,"m", "chord length of thruster")
e_d              = Variable("e_d", 0.202e-3,"m", "emitter diamaeter")
C_De             = Variable("C_De", 2,"-", "emitter drag coefficient")
A_t              = Variable("A_t", "m^2", "Thruster characteristic area")
L_t              = Variable("L_t", "m", "Thruster characteristic length")
Ke_fac           = Variable("Ke_fac", "-", "Multiplier for stage drag effects, emitter electrode")
Kg_fac           = Variable("Kg_fac", "-", "Multiplier for stage drag effects, ground electrode")
E_t              = Variable("E_t", "V/m", "Thruster electric field")
V_iw             = Variable("V_iw", "m/s", "ionic wind velocity")
G_t              = Variable("G_t", 0.05, "m", "Thruster gap distance")
T2P              = Variable("T2P", "N/kW", "Thruster thrust-to-power ratio")
V_op             = Variable("V_{op}", "kV", "Thruster operating voltage")
DeltaL           = Variable("DeltaL", "-", "non-dimensional pair spacing")
TA               = Variable("TA", "N/m^2", "thrust per unit area")
Trat             = Variable("Trat", "-", "thrust ratio")
T1               = Variable("T1", "N", "thrust of single pair")
Gain             = Variable("Gain", 200, "-", "Voltage Gain of transformer")
eta_prop         = Variable("eta_prop", "-", "Propulsive efficiency")

# Drag variables
C_Dp             = Variable("C_{D_p}", "-", "drag model parameter")
C_DA             = Variable("C_DA", "m^2", "Total vehicle drag times area")
C_D              = Variable("C_D", "-", "Total vehicle drag coefficient")
C_Dt             = Variable("C_Dt", "-", "Thruster drag coefficient")
C_Dss            = Variable("C_Dss", "-", "Side strut drag coefficient")
C_Di_w           = Variable("C_{D_iw}", "-", "induced drag coefficient, wing")
C_DA_F           = Variable("C_{D}A_{F}", "m^2", "Fuselage drag coefficient * frontal area")
C_D_ht           = Variable("C_{D_{ht}}", "-", "Horizontal Tail Drag coefficent")
C_Di_ht          = Variable("C_{D_{i_{HT}}}", "-", "induced drag coefficient, HT")
e                = Variable("e", 0.95, "-", "Wing spanwise efficiency")
Re               = Variable("Re", "-", "Reynold's number")
Re_f             = Variable("Re_{f}", "-", "Reynold's number, for fuselage (MTFLOW)")
Re_ss            = Variable("Re_{ss}", "-", "Reynold's number, for side strut")
Re_t             = Variable("Re_{t}", "-", "Reynold's number, for thruster")
D                = Variable("D", "N", "Total Drag on aircraft")
Dt               = Variable("Dt", "N", "Drag on thruster")
Dnt              = Variable("Dnt", "N", "Drag not on thruster")

# Structural variables
N_lift           = Variable("N_{lift}", "N", "Wing loading multiplier")
sig_maxh60       = Variable("\\sigma_{max,h60}", 15e3, "Pa", "Allowable stress, highload 60")
sigs_maxh60      = Variable("\\sigma_{max,shear,h60}", 400e3, "Pa", "Allowable shear stress")
sig_maxcf        = Variable("\\sigma_{max,cf}", 570e6, "Pa", "Allowable stress, carbon fibre")
sigs_maxcf       = Variable("\\sigma_{max,shear,cf}", 90e6, "Pa", "Allowable shear stress")
E_h60            = Variable("E_{h60}", 15000e3, "Pa", "Young's Modulus of high load 60 foam")
E_cf             = Variable("E_{cf}", 140e9, "Pa", "Young's Modulus of carbon fibre")
G_h60            = Variable("G_{h60}", 7.5e6, "Pa", "Shear Modulus high load 60 foam")
delta_max        = Variable("\\delta_{max}}", 0.05, "m", "Maximum tip deflection of wing")
Mom_w            = Variable("Mom_{w}", "Nm", "root bending moment of wing")
MomW_w           = Variable("MomW_{w}", "Nm", "root bending moment of wing, due to weight")
I_w              = Variable("I_{w}", "m^4", "second area moment of wing spar")
Q                = Variable("Q", "N", "shear force in wing spar")
QW               = Variable("QW", "N", "shear force in wing spar due to weight")
delw_b           = Variable("delw_{b}", "m", "deflection due to bending")
delw_s           = Variable("delw_{s}", "m", "deflection due to shear")
del_max          = Variable("del_{max}", 0.025, "-", "max % displacement allowed")
delw             = Variable("delw", "m", "total deflection of wing")

# Flight variables
V                = Variable("V", "m/s", "Flight speed")
C_L              = Variable("C_L", "-", "Wing lift coefficent")
L                = Variable("L", "N", "Total wing lift")
S                = Variable("S", "m^2", "Wing area")
A                = Variable("A", "-", "Aspect Ratio")
gam              = Variable("\\gamma","rad","Flight angle")
E                = Variable("E", "-", "Aerodynamic efficiency")

# Stability variables
xnp              = Variable("xnp","m","neutral point location")
xcg              = Variable("xcg","m","center of gravity location")
zcg              = Variable("zcg", "m", "center of gravity location in z-axis")
C_M              = Variable("C_{M}", -0.0440, "-", "Averaged Moment coefficient of tS12 airfoil")
C_M_ht           = Variable("C_{M_{HT}}", -0.0453, "-", "Averaged moment coefficient of HT")
htf              = Variable("htf", "-", "Vh*S/(l_ht*c)")

# Range equations variables
E_spec           = Variable("E_spec", "W*h/kg", "Mass specific energy content of battery")
Range            = Variable("Range", "m", "Range of UAV")
Endurance        = Variable("Endurance", "s", "Endurance of flight")
Rmin             = Variable("R_min", "m", "Minimum turn radius")

# Sweep variables
alfa1            = Variable("alfa_1", 0.4, "-", "Sweep coefficient")
alfa2            = Variable("alfa_2", 0.6667, "-", "Sweep coefficient")
alfa3            = Variable("alfa_3", 0.6667, "-", "Sweep coefficient")

# Utopian values for objective function
Range_ut         = Variable("Range_ut", 2459.09, "m", "Utopian Range")
E_ut             = Variable("E_ut", "-", 18.60, "Utopain aerodonyamic efficiency")
W_ut             = Variable("W_ut", 10.29, "N", "Utopian weight")
eta_prop_ut      = Variable("eta_prop_ut", 0.0505 ,"-", "Utopian propulsive efficiency")

gpkit.enable_signomials()

eqns     = [

#### Flight Model ####

                L == 0.5*rho*C_L*S*V**2*di_corr, # define lift
                Dnt == 0.5*rho*V**2*C_DA, # define drag, everything but thruster
                Dt >= 0.5*rho*V**2*(C_Dt*K*A_t+C_De*L_t*e_d) +
                          0.5*rho*V_iw**2*(C_Dt*Kg_fac*A_t+C_De*Ke_fac*L_t*e_d), # define drag, thruster
                D >= Dnt + Dt, # define total drag

                # Impose condition for Stall C_L vs Re relationship
                1 >= 176.8*C_L**17.5/Re + 7.429e5*C_L**2.275/Re - 7.425e5*C_L**2.129/Re,

                Re == (1/nu)*V*c, # define wing Re
                Re_f == (1/nu)*V*0.10668, # define fuselage Re (MTFLOW)
                Re_t == (1/nu)*V*2*0.0254, # define ground electrode Re
                Re_ss == (1/nu)*V*20*0.0254, # define side strut Re
                Re_ht == (1/nu)*V*c_ht, # define HT Re

                ## Steady Level Flight ##
                L >= W,
                T >= D,

                ## Climbing ##
#               T >= W*gam + D, 
#               L >= L,

                b == (S*A)**0.5, # define wing span
                c == (S/A)**0.5, # define wing chord
                C_D == D/(0.5*rho*V**2*S), # calculate vehicle drag coefficient  

                E == L/D, # Aerodynamic efficiency

#### Thrust ####

                T == (A*S)**0.5*ndi_frac*N*K/(T1*Trat), # get thrust given thruster geometry
                T1 >= 1221695966.01*V_op**-7.3747673794*1.e3/5.,
                T1 >= 4151.46592088*V_op**-3.12069777202*1.e3/5.,
                T1 >= 67112.3207544*V_op**-3.98906732575*1.e3/5.,
                Trat**54.74862 >= 0.06725*DeltaL**-13.99422 + 0.91309*DeltaL**-1.92270,
                P == T/T2P, # estimate power from T/P
                P <= 0.750, # limited by expected output from PC prototype
                E_t == V_op*1.e3/G_t, # calculate thruster electric field (1D)
                V_iw == (E_t**2*eps0/rho)**0.5, # calculate ionic wind velocity
                K >= 1,
                N >= 1,
                K <= 3,
                T2P == G_t/(2e-4*V_op),
                DeltaL*N*G_t <= t_f,
                V_op >= 15,
                V_op <= 40,
                n_batt == V_op*1000/(Gain*3.3),
                DeltaL <= 2.0,
                DeltaL >= 0.4,

                eta_prop == T*V/(P*1000),

#### Drag Model ####

                C_DA_F**2 == 4.897e-2*Re_f**-0.9811, # fuselage drag

                C_Di_w >= C_L**2/(pi*A*e) + 0.38*C_Dp*C_L**2, # wing induced drag (Schaufelle model)

                C_Di_ht >= C_L_ht**2/(pi*A*e) + 0.38*C_D_ht*C_L_ht**2, # HT induced drag (Schaufelle model)

                C_Dss == 1.33/Re_ss**0.5, # side strut drag

                1 >= (0.0309*C_L**4.9014/(Re**1.9359*C_Dp**5.826) +
                     166.9188*C_L**20.4236/Re**3.4473*C_Dp**8.1109) +
                     0.1722*Re**0.0887/(C_L**0.0722*C_Dp**0.0957) +
                     1.8286e7*C_L**0.1999/(Re**2.4252*C_Dp**2.3477), # wing profile drag

                A_t == ndi_frac*(A*S)**0.5*c_t*N, # planform area, ground electrodes
                L_t == ndi_frac*(A*S)**0.5*N, # length of electrodes per stage
                Kg_fac == 0.921594*K**1.75391, # factor to account for increasing speed between stages
                Ke_fac == Kg_fac, # fit no good, assume this as "worst case"

                # thruster drag
                C_Dt**10.5904 >= 5.6706e-9*Re_t**-1.5691 + 4750.3152*Re**-4.4195 + 7.4263e-6*Re**-2.9821,

                # vehicle drag coeff. times area, no thruster
                C_DA >= C_DA_F  + C_Dp*S + C_Di_w*S + 2*2*C_Dss*A_ss + C_Di_ht*S_ht             + C_D_ht*S_ht,

#### Weight Constraints ####

                W >=  W_f + W_w + W_batt + W_pc + W_t + W_ht + W_vt + W_ts + W_ss, # total weight
                P/eta_pc/n_batt <= 15.e-3, # limit cell power
                W_batt == g*n_batt*Mb, # battery weight
                nr == b/rs, # calculate number of ribs in wing
                W_ml == g*SA_ref*c*b*ml_sd, # weight of microlite
                W_le == g*Ale_ref*c**2*b*rho_f250, # weight of LE
                W_te == g*Ate_ref*c**2*b*rho_f250, # weight of TE
                W_r == g*Ar_ref*c**2*rt*nr*rho_f250, # weight of ribs
                W_s >= g*(c**2*tau*st)*b*rho_h60 + g*cft*(st*c)*b*rho_cf, # weight of spar
                W_w >= W_ml + W_te + W_le + W_r + W_s, # wing weight
                W_t == g*0.0134*(A*S)**0.5*ndi_frac*N*K, # thruster weight
                W_ts == g*l_ht*tspm, # tail spar weight
                W_ht == W_w*S_ht/S, # HT weight
                W_vt == W_w*S_vt/S, # VT weight
                W_pc == g*0.375, # power converter weight
                W_ss == 2*g*sspm2*A_ss, # side strut weight

#### Structural Model ####

                Mom_w == L/2*b/4, # calculate root bending moment, from lift (worst case, in flight)
                # calculate root bending moment, from weight (not in flight, worst case)
                MomW_w >= W_w/2*b/4 + W_ss/2*ndi_frac*b/2 + W_t/2*ndi_frac*b/4,
                I_w == 2*(st*c)*cft*(tau*c/2)**2, # spar moment of inertia (just cf)
                sig_maxcf >= (N_lift*Mom_w*c*tau/2)/I_w, # bending stress limit
                sig_maxcf >= (N_lift*MomW_w*c*tau/2)/I_w, # bending stress limit, weight
                Q == L/2, # calculate root shear stress
                QW >= W_w/2 + W_ss/2 + W_t/2, # calculate root shear stress due to weight
                sigs_maxh60 >= N_lift*Q/(st*tau*c**2), # shear stress limit (just foam)
                sigs_maxh60 >= N_lift*QW/(st*tau*c**2), # shear stress limit (just foam), weight
                delw_b == N_lift*L/2*(b/2)**3/(8*E_cf*I_w), # deflection due to bending
                delw_s == N_lift*Mom_w/(2*G_h60*(st*tau*c**2)), # deflection due to shear
                delw >= delw_b + delw_s, # total deflection
                delw/b <= del_max, # deflection limit

#### Turning requirements ####

                N_lift == 0.5*rho*V**2*C_L_max/(W/S), # C_L_max for tS12 = 1.5312

                Rmin == 1.45*V**2/(g*N_lift**1.299), # minimum turn radius (posynomial approximation)

                Rmin <= 50/pi, # minimum turning radius must be smaller than track field radius

#### Stability requirements ####

# get required tail lift coefficient

                # Moment arms for now
                zcg == 0.01, # must be defined (cannot = 0)
                l_w == 0.1,

                # Lift coefficient requirement
                C_L >= C_L_ht,

                # xcp_w < xcg
                0 >= l_w*C_L/(C_L_ht*htf*l_ht) + c*C_M/(C_L_ht*htf*l_ht) + c_ht*C_M_ht*htf/(C_L_ht*htf*l_ht) + zcg*T/L/(C_L_ht*htf*l_ht) - 1,

# get drag coefficient of HT

                1 >= 0.2201*C_D_ht**-0.2762*C_L_ht**-0.078*Re_ht**-0.3878*tau_ht**0.8912
                        + 0.0046*C_D_ht**-0.2109*C_L_ht**0.472*Re_ht**0.2637*tau_ht**-0.1398
                        + 2.1243*C_D_ht**-0.4001*C_L_ht**-0.0241*Re_ht**-0.2605*tau_ht**0.1337,

                tau_ht >= 0.08,
                tau_ht <= 0.14,

#### Tail properties ####

# get HT span, chord, and area
                S_ht == Vh*S*c/l_ht,
                b_ht == (A_ht*S_ht)**0.5,
                c_ht == S_ht/b_ht,

                htf == Vh*c/(l_ht),

# get VT span, chord, and area
                S_vt == Vv*S*b/l_vt,
                b_vt == (A_vt*S_vt)**0.5,
                c_vt == S_vt/b_vt,
#### COG Calculations ####

# get neutral point

                (c/xnp) == 2.18270636553*A**-0.106088939277,

#### Range Equations ####

                E_spec == (0.33)/Mb, # Wh/kg
                Range == V*(W_batt/g)*E_spec/(P*1000)*60**2,
                Endurance == Range/V,

                ]

#### Objective function ####

Obj = alfa1*alfa2*alfa3*(W/W_ut) + alfa2*alfa3*(E_ut/E) - alfa1*alfa2*alfa3*(E_ut/E) + alfa3*(Range_ut/(Range)) - alfa2*alfa3*(Range_ut/(Range)) + (eta_prop_ut / (T*V/(P*1000))) - alfa3*(eta_prop_ut / (T*V/(P*1000)))

m = gpkit.Model(Obj, eqns)
sol = m.localsolve()
whoburg commented 8 years ago

Thanks @bs1912 ! As @bqpd mentioned above, we're going to keep this open so we can use this example to motivate some internal improvements. Thanks again.

bqpd commented 8 years ago

Alright, so one way to solve this problem (sig constraints only solve when of a particular form, which used to be default and no longer is) would be to turn signomial constraints of the form m >= s into 0 >= s/m -1, as we used to. It's kind of always a bad idea to take your biggest term and make it negative (and thus requiring approximation). Another would be to mention this in the docs as a hazard of signomials.

whoburg commented 8 years ago

Is there any chance we're going to accidentally resolve this by developing better SP heuristics? I'm not clear on the exact reason the new heuristic fails on the 0 >= s - m interpretation.

bqpd commented 8 years ago

Yup, this is all about better SP heuristics. Specifically the way in which our previous one was better in this case, because it kept the largest term in the posynomial (which as we've discussed before is helpful for a good negynomial approximation).

whoburg commented 8 years ago

Which specific constraint had the change in approximation?

bqpd commented 8 years ago

All the ones that start with 1 >= above

whoburg commented 8 years ago

It's just the Reynolds number constraint, right?

In the new way, that constraint is Re + 7.425e5/C_L**-2.129 >= 176.8/C_L**-17.5 + 7.429e5/C_L**-2.275.

Does our current heuristic set the initial guess for both Re and CL to 1? That's a pretty awful guess for Re...

For this case I think the right way is to "treat signomials as posynomial constraints, ignoring negative terms", as @bqpd said. I suspect that's different (and better!) than going back to the old 1 >= way.

bqpd commented 8 years ago

Hmm, but that old approximation method couldn't solve all problems either...

whoburg commented 8 years ago

True. Do we agree the fundamental failure here has to do with the [1, 1] initial guess?

bqpd commented 8 years ago

That plus that fact that 0 >= changes the significance of that guess...there's still two parts to this, though I am on board with 0 >= from a mathematical perspective, it's caused us mostly trouble from a practical one.

whoburg commented 8 years ago

I also wonder whether the "all monomial terms in a negynomial have equal contributions" heuristic for choosing an initial guess would work here.

bqpd commented 7 years ago

I think I'm going to close this one and opens a Docs issue noting that SP formulation does affect solving.