generating united constants for power laws

[codykarcher]

So for weights, I have lots of empirical trends that look something like this:

S          = Variable("S", 800, "ft^2", "Wing Area")
k_w        = Variable("k_w", 1.7e-3, "-", "Proportionality factor for wing")
n_ult      = Variable("n_{ult}", 2.5, "-", "Maximum Load Factor")
MTOW  = Variable("W_{MTO}", "lbf", "Maximum Takeoff Weight")
bs_ratio   = Variable("bs_ratio", 1.2, "-", "Ratio of unswept span to span")
b_ref      = Variable("b_{ref}", 200, "ft", "Wing Span")
t_r        = Variable("\\tau_{root}", .15, "-", "Root Chord Thickness")

W_wing >= MTOW * k_w * b_ref**(0.75) * (1+(bs_ratio)**(0.5)) * n_ult**(0.55) * ((b_ref/t_r)/(MTOW/S))**(0.30)

which is a spectacular units nightmare, however, this is the trend that is used. What I want to be able to do is disable the units only for this constraint, as opposed to having to back out what the units are and then add crazy unit terms in order to cancel it out and make it work. Is this possible?

Note: Don't use this example for anything practical. I messed with some stuff to make it postable.

[codykarcher]

What I really want is for the gpkit.disable_units() environment to have similar functionality to the gpkit.Signomials_Enabled environment.

[bqpd]

Haha. We certainly could do that. @whoburg, any objections?

[whoburg]

nope, this sounds reasonable (and I agree with @cjk12 that it's common to find power law relationships given where it would be a pain to figure out the units of the constants).

Couple items to make sure we're considering: Are there any issues expected with using a constraint constructed in this way in other contexts? E.g. is any unit checking done post-constraint creation? And are there any places where we rely on the fact that things are guaranteed to have the right units?

[codykarcher]

What's the priority on this? I could really use it now that I've got the Raymer weights model working.

[whoburg]

I'd budget for this talking a little while - we need to make sure we don't screw it up.

For example, if I make a constraint set S involving one of these constraints with disabled units, and then use it in a large model that has different units on one or more variables in S, is there a possibility there will be an uncaught units mismatch?

@cjk12, can I see an example of a power law you want this for? I'm wondering whether some of these can be made unitless by being careful to make each term non-dimensional.

[codykarcher]

Here is my weights breakdown:

# Bradley
W_centerbody >= 5.698865 * 0.316422 * (W_maximum_gross ** 0.166552) * Scab ** 1.061158,
W_afterbody >= (1.0 + 0.05*Nen) * 0.53 * Saft * (W_maximum_gross**0.2) * (lambda_aft + 0.5),
# Raymer page 459
W_wing >= W_maximum_gross * Kw * (Bw)**(0.75) * (1+(Bref/Bw)**(0.5)) * Nz**(0.55) * (Bw/(tau_w*chord))**(0.3) * (W_maximum_gross/(Sw*.2))**(-0.3), #Torenbeek for GP compatibility
W_vertical_tail >= 12000,
W_main_landing_gear >= 0.0106 * Kmp * W_landing_gross**(0.888) * Nl**(0.25) * Lm**(0.4) * Nmw**(0.321) * Nmss**(-0.5) * V_stall**(0.1),
W_nose_landing_gear >= 0.032 * Knp * W_landing_gross**(0.646) * Nl**(0.2) * Ln**(0.5) * Nnw**(0.45),
W_engine_contents >= 2.331 * Wen**(0.901) * Kp * Ktr,
W_nacelle_group >= 0.6724 * Kng * NLt**(0.10) * Nw**(0.294) * Nz**(0.119) * W_engine_contents**(0.611) * Nen**(0.984) * Sn**(0.224),
W_engine_controls >= 5.0 * Nen + 0.80 * Lec,
W_starter_penumatic >= 49.19 * (Nen * Wen / 1000.)**(0.541),
W_fuel_system >= 2.405 * (Vt*7.48052)**(0.606) * Nt**(0.5),
W_flight_controls >= 0.005 * W_maximum_gross,
W_APU_installed >= 2.2 * W_APU_uninstalled,
W_instruments >= Kieg * W_operating_empty**(0.555) * R**(0.25), # Torenbeek for GP compatibility
W_hydraulics >= 0.2673 * Nf * (Lf + Bw),
W_electrical >= 7.291 * Rkva**(0.782) * La**(0.346) * Ngen**(0.10),
W_avionics >= 1.73 * Wuav**(0.983),
W_furnishings >= 0.0577 * Nc**(0.1) * Wc**(0.393) * Scab**(0.75),
W_air_conditioning >= 62.36 * Np**(0.25) * (Vpr/1000.)**(0.604) * Wuav**(0.01),
W_handling_gear >= 0.002 * W_maximum_gross,
W_anti_ice >= 0.0003 * W_maximum_gross,
W_military_cargo_handling_system  >= 2.4 * Scab,

I can definitely do the units by hand, (I've done it before), but was just hoping to avoid it, since it's a lot of slow debugging (see #657). I'll probably just set up units constants or something which cancel out the units.

[whoburg]

Yeah this is a very real problem but let's think carefully about how to solve it...

The power laws listed in these textbooks are only valid for certain units. At the very least, the code should say exactly what the right units are (and perhaps assert that the variables have the right units?)

An even better way would be to try to write these in non-dimensional form, as Drela suggests. For example, you would aim for something like

W_centerbody/W_cb_ref >= (some unitless constant) * (W_maximum_gross/W_maxgross_ref)**some_exp1 * (Scab/Scab_ref)**some_exp2.

This way, the reference quantities encode the units the power law was fitted for -- so W_maxgross_ref can be in lbs and W_maximum_gross in newtons, and all will work fine.

I wonder whether instead of disabling units for a particular constraint, we should write a converter that converts power laws in the form above (with correct units given as input) to the dimensionless reference quantity form. Seems way more robust.

[codykarcher]

Yea, that's definitely an option. Couple complicating factors, first, most of the leading constants aren't actually unitless, and switch depending on what set of units are being used. I can probably work it out though. Second, there are also constraints like: W_engine_controls >= 5.0 * Nen + 0.80 * Lec which add things that don't have the same units... This can also obviously be fixed by hand, but may cause problems for an automatic routine.

[whoburg]

What are the units of W_engine_controls, Nen , and Lec?

[codykarcher]

pounds, unitless and feet respectively

[whoburg]

cool, so one way to write the above is

W_engine_controls/(5 * gpkit.units("lbf")) >= Nen + Lec/(6.25 * gpkit.units("ft")). That constraint is robust to whatever units of W_engine_controls or Lec you throw at it.

It's of course not unique -- I could have instead written W_engine_controls/gpkit.units("lbf") >= 5.0*Nen + Lec/(1.25 * gpkit.units("ft")).

I claim that this is always possible to do, even in power laws with ugly exponents.

If the power law pertains to a certain size vehicle or subsystem, a good way to resolve the ambiguity of what to divide by is to make the reference quantity for that parameter equal to a representative value -- e.g. if the model above tends to be for Lec around 6 feet, then make Lec_ref equal 6 feet. If it's valid when W_engine_controls is near 10 lb, then make W_engine_controls_ref equal 10lb.

[whoburg]

I also think this can be automated. The inputs are

1. the textbook power law
2. the units of all quantities
3. one or more values chosen as the reference values, depending on the multiplicity of the monomial terms

[bqpd]

@whoburg, would this be something for the tools submodule?

[bqpd]

@whoburg, I think this is pretty low priority?

[whoburg]

It's low priority, but relevant enough to GP modeling that it should definitely stay open.

[pgkirsch]

Instead of generating unit-ed constants automatically would it be possible to introduce a .unitless class method as a simple solution? My current way of handling the problem of power law/fitted models is (similar to @whoburg's suggestion on April 17th):

z >= 3.56*z.units*(x/x.units)**0.983*(y/y.units)**1.234,

This can make for very long constraints if x and y are more descriptive variable names. Suggested solution (and something that would be nice to implement anyway):

z >=3.56*z.units*(x.unitless)**0.983*(y.unitless)**1.234,

[bqpd]

short answer: EWW

long answer:

z >= 3.56*z.units*(x/gpkit.units("m"))**0.983*(y/gpkit.units("volts"))**1.234,

I'd be willing to write a .unitless_in_SI_base_units method, but not a generic unitless method.

[pgkirsch]

Yeah sure, but that's (very slightly) less robust/general. Any reason a unitless method would make you sad?

[pgkirsch]

Sorry I'm not sure I understand why unitless_in_SI_base_units makes you happier.. what am I missing?

[bqpd]

On the contrary, since it gets the same answer even if x and y are linked with a variable with different-but-compatible units, I'd say it's much more robust.

(linking being why these unitless methods make me sad)

[whoburg]

Yeah, z >= 3.56*z.units*(x/x.units)**0.983*(y/y.units)**1.234 scares me -- that's only valid for some specific units of x, y, and z, but will compile fine independent of units. I'd agree that it's better to hard-code the units that the power law assumes.

[bqpd]

Discussion was in favour of harcoded units.