Building Complex Models

Inheriting from Model

GPkit encourages an object-oriented modeling approach, where the modeler creates objects that inherit from Model to break large systems down into subsystems and analysis domains. The benefits of this approach include modularity, reusability, and the ability to more closely follow mental models of system hierarchy. For example: two different models for a simple beam, designed by different modelers, should be able to be used interchangeably inside another subsystem (such as an aircraft wing) without either modeler having to write specifically with that use in mind.

When you create a class that inherits from Model, write a .setup() method to create the model’s variables and return its constraints. GPkit.Model.__init__ will call that method and automatically add your model’s name and unique ID to any created variables.

Variables created in a setup method are added to the model even if they are not present in any constraints. This allows for simplistic ‘template’ models, which assume constant values for parameters and can grow incrementally in complexity as those variables are freed.

At the end of this page a detailed example shows this technique in practice.

Accessing Variables in Models

GPkit provides several ways to access a Variable in a Model (or ConstraintSet):

  • using Model.variables_byname(key). This returns all Variables in the Model, as well as in any submodels, that match the key.
  • using Model.topvar(key). This returns the top-level Variable that matches the key. The Variable must appear at the top level, not in a submodel.
  • using Model.__getitem__. Model[key] returns the only variable matching the key, even if the match occurs in a submodel. If multiple variables match the key, an error is raised.

These methods are illustrated in the following example.

"Demo of accessing variables in models"
from gpkit import Model, Variable


class Battery(Model):
    "A simple battery"
    def setup(self):
        h = Variable("h", 200, "Wh/kg", "specific energy")
        E = Variable("E", "MJ", "stored energy")
        m = Variable("m", "lb", "battery mass")
        return [E <= m*h]


class Motor(Model):
    "Electric motor"
    def setup(self):
        m = Variable("m", "lb", "motor mass")
        f = Variable("f", 20, "lb/hp", "mass per unit power")
        Pmax = Variable("P_{max}", "hp", "max output power")
        return [m >= f*Pmax]


class PowerSystem(Model):
    "A battery powering a motor"
    def setup(self):
        components = [Battery(), Motor()]
        m = Variable("m", "lb", "mass")
        return [components,
                m >= sum(comp.topvar("m") for comp in components)]

PS = PowerSystem()
print "Getting the only var 'E': ", PS["E"]
print "The top-level var 'm': ", PS.topvar("m")
print "All the variables 'm': ", PS.variables_byname("m")
Getting the only var 'E':  E_PowerSystem/Battery [MJ]
The top-level var 'm':  m_PowerSystem [lb]
All the variables 'm':  [gpkit.Variable(m_PowerSystem [lb]), gpkit.Variable(m_PowerSystem/Battery [lb]), gpkit.Variable(m_PowerSystem/Motor [lb])]

Vectorization

gpkit.Vectorize creates an environment in which Variables are created with an additional dimension:

"from gpkit/tests/t_vars.py"

def test_shapes(self):
    with gpkit.Vectorize(3):
        with gpkit.Vectorize(5):
            y = gpkit.Variable("y")
            x = gpkit.VectorVariable(2, "x")
        z = gpkit.VectorVariable(7, "z")

    self.assertEqual(y.shape, (5, 3))
    self.assertEqual(x.shape, (2, 5, 3))
    self.assertEqual(z.shape, (7, 3))

This allows models written with scalar constraints to be created with vector constraints:

"Vectorization demonstration"
from gpkit import Model, Variable, Vectorize

class Test(Model):
    "A simple scalar model"
    def setup(self):
        x = Variable("x")
        return [x >= 1]

print "SCALAR"
m = Test()
m.cost = m["x"]
print m.solve(verbosity=0).summary()

print "__________\n"
print "VECTORIZED"
with Vectorize(3):
    m = Test()
m.cost = m["x"].prod()
m.append(m["x"][1] >= 2)
print m.solve(verbosity=0).summary()
SCALAR

Cost
----
 1

Free Variables
--------------
x : 1

__________

VECTORIZED

Cost
----
 2

Free Variables
--------------
x : [ 1         2         1        ]

Multipoint analysis modeling

In many engineering models, there is a physical object that is operated in multiple conditions. Some variables correspond to the design of the object (size, weight, construction) while others are vectorized over the different conditions (speed, temperature, altitude). By combining named models and vectorization we can create intuitive representations of these systems while maintaining modularity and interoperability.

In the example below, the models Aircraft and Wing have a .dynamic() method which creates instances of AircraftPerformance and WingAero, respectively. The Aircraft and Wing models create variables, such as size and weight without fuel, that represent a physical object. The dynamic models create properties that change based on the flight conditions, such as drag and fuel weight.

This means that when an aircraft is being optimized for a mission, you can create the aircraft (AC in this example) and then pass it to a Mission model which can create vectorized aircraft performance models for each flight segment and/or flight condition.

"""Modular aircraft concept"""
import numpy as np
from gpkit import Model, Variable, Vectorize


class Aircraft(Model):
    "The vehicle model"
    def setup(self):
        self.fuse = Fuselage()
        self.wing = Wing()
        self.components = [self.fuse, self.wing]

        W = Variable("W", "lbf", "weight")

        return self.components, [
            W >= sum(c.topvar("W") for c in self.components)
            ]

    def dynamic(self, state):
        "This component's performance model for a given state."
        return AircraftP(self, state)


class AircraftP(Model):
    "Aircraft flight physics: weight <= lift, fuel burn"
    def setup(self, aircraft, state):
        self.aircraft = aircraft
        self.wing_aero = aircraft.wing.dynamic(state)
        self.perf_models = [self.wing_aero]
        Wfuel = Variable("W_{fuel}", "lbf", "fuel weight")
        Wburn = Variable("W_{burn}", "lbf", "segment fuel burn")

        return self.perf_models, [
            aircraft.topvar("W") + Wfuel <= (0.5*state["\\rho"]*state["V"]**2
                                             * self.wing_aero["C_L"]
                                             * aircraft.wing["S"]),
            Wburn >= 0.1*self.wing_aero["D"]
            ]


class FlightState(Model):
    "Context for evaluating flight physics"
    def setup(self):
        Variable("V", 40, "knots", "true airspeed")
        Variable("\\mu", 1.628e-5, "N*s/m^2", "dynamic viscosity")
        Variable("\\rho", 0.74, "kg/m^3", "air density")


class FlightSegment(Model):
    "Combines a context (flight state) and a component (the aircraft)"
    def setup(self, aircraft):
        self.flightstate = FlightState()
        self.aircraftp = aircraft.dynamic(self.flightstate)
        return self.flightstate, self.aircraftp


class Mission(Model):
    "A sequence of flight segments"
    def setup(self, aircraft):
        with Vectorize(4):  # four flight segments
            self.fs = FlightSegment(aircraft)

        Wburn = self.fs.aircraftp["W_{burn}"]
        Wfuel = self.fs.aircraftp["W_{fuel}"]
        self.takeoff_fuel = Wfuel[0]

        return self.fs, [Wfuel[:-1] >= Wfuel[1:] + Wburn[:-1],
                         Wfuel[-1] >= Wburn[-1]]


class Wing(Model):
    "Aircraft wing model"
    def dynamic(self, state):
        "Returns this component's performance model for a given state."
        return WingAero(self, state)

    def setup(self):
        W = Variable("W", "lbf", "weight")
        S = Variable("S", 190, "ft^2", "surface area")
        rho = Variable("\\rho", 1, "lbf/ft^2", "areal density")
        A = Variable("A", 27, "-", "aspect ratio")
        c = Variable("c", "ft", "mean chord")

        return [W >= S*rho,
                c == (S/A)**0.5]


class WingAero(Model):
    "Wing aerodynamics"
    def setup(self, wing, state):
        CD = Variable("C_D", "-", "drag coefficient")
        CL = Variable("C_L", "-", "lift coefficient")
        e = Variable("e", 0.9, "-", "Oswald efficiency")
        Re = Variable("Re", "-", "Reynold's number")
        D = Variable("D", "lbf", "drag force")

        return [
            CD >= (0.074/Re**0.2 + CL**2/np.pi/wing["A"]/e),
            Re == state["\\rho"]*state["V"]*wing["c"]/state["\\mu"],
            D >= 0.5*state["\\rho"]*state["V"]**2*CD*wing["S"],
            ]


class Fuselage(Model):
    "The thing that carries the fuel, engine, and payload"
    def setup(self):
        # fuselage needs an external dynamic drag model,
        # left as an exercise for the reader
        # V = Variable("V", 16, "gal", "volume")
        # d = Variable("d", 12, "in", "diameter")
        # S = Variable("S", "ft^2", "wetted area")
        # cd = Variable("c_d", .0047, "-", "drag coefficient")
        # CDA = Variable("CDA", "ft^2", "drag area")
        Variable("W", 100, "lbf", "weight")

AC = Aircraft()
MISSION = Mission(AC)
M = Model(MISSION.takeoff_fuel, [MISSION, AC])
sol = M.solve(verbosity=0)

vars_of_interest = set(AC.varkeys)
vars_of_interest.update(MISSION.fs.aircraftp.unique_varkeys)
vars_of_interest.add("D")
print sol.summary(vars_of_interest)

Note that the output table can be filtered with a list of variables to show.


Cost
----
 1.943 [lbf] 

Free Variables
--------------
         | Aircraft
       W : 290                                         [lbf] weight

         | Aircraft/Wing
       W : 190                                         [lbf] weight
       c : 2.653                                       [ft]  mean chord

         | Mission/FlightSegment/AircraftP
W_{burn} : [ 0.487     0.486     0.485     0.485    ]  [lbf] segment fuel burn
W_{fuel} : [ 1.94      1.46      0.97      0.485    ]  [lbf] fuel weight

         | Mission/FlightSegment/AircraftP/WingAero
       D : [ 4.87      4.86      4.85      4.85     ]  [lbf] drag force

Sensitivities
-------------
     | Aircraft/Fuselage
   W : +0.25            weight

     | Aircraft/Wing
   S : +0.68            surface area
\rho : +0.48            areal density
   A : -0.31            aspect ratio

Next Largest Sensitivities
--------------------------
     | Mission/FlightSegment/AircraftP/WingAero
   e : [ -0.093    -0.092    -0.092    -0.092   ] Oswald efficiency

     | Mission/FlightSegment/FlightState
   V : [ +0.1      +0.1      +0.1      +0.1     ] true airspeed
\rho : [ +0.034    +0.034    +0.035    +0.035   ] air density