3

I'm interested in using SymPy to augment my engineering models. Instead of defining a rigid set of inputs and outputs, I'd like for the user to simply provide everything they know about a system, then apply that data to an algebraic model which calculates the unknowns (if it has enough data).

For an example, say I have some stuff which has some mass, volume, and density. I'd like to define a relationship between those parameters (density = mass / volume) such that when the user has provided enough information (any 2 variables), the 3rd variable is automatically calculated. Finally, if any value is later updated, one of the other values should change to preserve the relationship. One of the challenges with this system is that when there are multiple independent variables, there would need to be a way to specify which independent variable should change to satisfy the requirement.

Here's some working code that I have currently:

from sympy import *

class Stuff(object):
    def __init__(self, *args, **kwargs):
        #String of variables
        varString = 'm v rho'

        #Initialize Symbolic variables
        m, v, rho = symbols(varString)

        #Define density equation
        # "rho = m / v" becomes "rho - m / v = 0" which means the eqn. is rho - m / v
        # This is because solve() assumes equation = 0
        density_eq = rho - m / v

        #Store the equation and the variable string for later use
        self._eqn = density_eq
        self._varString = varString

        #Get a list of variable names
        variables = varString.split()

        #Initialize parameters dictionary
        self._params = {name: None for name in variables}

    @property
    def mass(self):
        return self._params['m']

    @mass.setter
    def mass(self, value):
        param = 'm'
        self._params[param] = value
        self.balance(param)

    @property
    def volume(self):
        return self._params['v']

    @volume.setter
    def volume(self, value):
        param = 'v'
        self._params[param] = value
        self.balance(param)

    @property
    def density(self):
        return self._params['rho']

    @density.setter
    def density(self, value):
        param = 'rho'
        self._params[param] = value
        self.balance(param)

    def balance(self, param):

        #Get the list of all variable names
        variables = self._varString.split()

        #Get a copy of the list except for the recently changed parameter
        others = [name for name in variables if name != param]

        #Loop through the less recently changed variables
        for name in others:
            try:
                #Symbolically solve for the current variable
                eq = solve(self._eqn, [symbols(name)])[0]
                #Get a dictionary of variables and values to substitute for a numerical solution (will contain None for unset values)
                indvars = {symbols(n): self._params[n] for n in variables if n is not name}
                #Set the parameter with the new numeric solution
                self._params[name] = eq.evalf(subs=indvars)
            except Exception, e:
                pass


if __name__ == "__main__":

    #Run some examples - need to turn these into actual tests

    stuff = Stuff()
    stuff.density = 0.1
    stuff.mass = 10.0
    print stuff.volume

    stuff = Water()
    stuff.mass = 10.0
    stuff.volume = 100.0
    print stuff.density

    stuff = Water()
    stuff.density = 0.1
    stuff.volume = 100.0
    print stuff.mass

    #increase volume
    print "Setting Volume to 200"
    stuff.volume = 200.0
    print "Mass changes"
    print "Mass {0}".format(stuff.mass)
    print "Density {0}".format(stuff.density)

    #Setting Mass to 15
    print "Setting Mass to 15.0"
    stuff.mass = 15.0
    print "Volume changes"
    print "Volume {0}".format(stuff.volume)
    print "Density {0}".format(stuff.density)

    #Setting Density to 0.5
    print "Setting Density to 0.5"
    stuff.density = 0.5
    print "Mass changes"
    print "Mass {0}".format(stuff.mass)
    print "Volume {0}".format(stuff.volume)

    print "It is impossible to let mass and volume drive density with this setup since either mass or volume will " \
          "always be recalculated first."

I tried to be as elegant as I could be with the overall approach and layout of the class, but I can't help but wonder if I'm going about it wrong - if I'm using SymPy in the wrong way to accomplish this task. I'm interested in developing a complex aerospace vehicle model with dozens/hundreds of interrelated properties. I'd like to find an elegant and extensible way to use SymPy to govern property relationships across the vehicle before I scale up from this fairly simple example.

I'm also concerned about how/when to rebalance the equations. I'm familiar with PyQt Signals and Slots which is my first thought for how to link dependent things together to trigger a model update (emit a signal when a value updates, which will be received by rebalancing functions for each system of equations that relies on that parameter?). Yeah, I really don't know the best way to do this with SymPy. Might need a bigger example to explore systems of equations.

Here's some thoughts on where I'm headed with this project. Just using mass as an example, I would like to define the entire vehicle mass as the sum of the subsystem masses, and all subsystem masses as the sum of component masses. In addition, mass relationships will exist between certain subsystems and components. These relationships will drive the model until more concrete data is provided. So if the default ratio of fuel mass to total vehicle mass is 50%, then specifying 100lb of fuel will size the vehicle at 200lb. However if I later specify that the vehicle is actually 210lb, I would want the relationship recalculated (let it become dependent since fuel mass and vehicle mass were set more recently, or because I specified that they are independent variables or locked or something). The next problem is iteration. When circular or conflicting relationships exist in a model, the model must be iterated on to hopefully converge on a solution. This is often the case with the vehicle mass model described above. If the vehicle gets heavier, there needs to be more fuel to meet a requirement, which causes the vehicle to get even heavier, etc. I'm not sure how to leverage SymPy in these situations.

Any suggestions?

PS Good explanation of the challenges associated with space launch vehicle design.

Edit: Changing the code structure based on goncalopp's suggestions...

class Balanced(object):
    def __init__(self, variables, equationStr):
        self._variables = variables
        self._equation = sympify(equationStr)

        #Initialize parameters dictionary
        self._params = {name: None for name in self._variables}

        def var_getter(varname, self):
            return self._params[varname]

        def var_setter(varname, self, value):
            self._params[varname] = value
            self.balance(varname)

        for varname in self._variables:
            setattr(Balanced, varname, property(fget=partial(var_getter, varname),
                                                fset=partial(var_setter, varname)))

    def balance(self, recentlyChanged):

        #Get a copy of the list except for the recently changed parameter
        others = [name for name in self._variables if name != recentlyChanged]

        #Loop through the less recently changed variables
        for name in others:
            try:
                eq = solve(self._equation, [symbols(name)])[0]
                indvars = {symbols(n): self._params[n] for n in self._variables if n != name}
                self._params[name] = eq.evalf(subs=indvars)
            except Exception, e:
                pass


class HasMass(Balanced):
    def __init__(self):
        super(HasMass, self).__init__(variables=['mass', 'volume', 'density'],
                                      equationStr='density - mass / volume')


class Prop(HasMass):
    def __init__(self):
        super(Prop, self).__init__()

if __name__ == "__main__":

    prop = Prop()
    prop.density = 0.1
    prop.mass = 10.0
    print prop.volume

    prop = Prop()
    prop.mass = 10.0
    prop.volume = 100.0
    print prop.density

    prop = Prop()
    prop.density = 0.1
    prop.volume = 100.0
    print prop.mass

What this makes me want to do is use multiple inheritance or decorators to automatically assign physical inter-related properties to things. So I could have another class called "Cylindrical" which defines radius, diameter, and length properties, then I could have like class DowelRod(HasMass, Cylindrical). The really tricky part here is that I would want to define the cylinder volume (volume = length * pi * radius^2) and let that volume interact with the volume defined in the mass balance equations... such that, perhaps mass would respond to a change in length, etc. Not only is the multiple inheritance tricky, but automatically combining relationships will be even worse. This will get tricky very quickly. I don't know how to handle systems of equations yet, and it's clear with lots of parametric relationships that parameter locking or specifying independent/dependent variables will be necessary.

Community
  • 1
  • 1
flutefreak7
  • 2,321
  • 5
  • 29
  • 39
  • Unfortunately, you can't assume solve will return a list. You should use `solve(dict=True)` to get a more consistent output. – asmeurer Mar 08 '14 at 00:59

3 Answers3

1

While I have no experience with this kind of models, and little experience with SymPy, here's some tips:

@property
def mass(self):
    return self._params['m']

@mass.setter
def mass(self, value):
    param = 'm'
    self._params[param] = value
    self.balance(param)

@property
def volume(self):
    return self._params['v']

@volume.setter
def volume(self, value):
    param = 'v'
    self._params[param] = value
    self.balance(param)

As you've noticed, you're repeating a lot of code for each variable. This is unnecessary, and since you'll have lots of variables, this eventually leads to code maintenance nightmares. You have your variables neatly arranged on varString = 'm v rho' My suggestion is to go further, and define a dictionary:

my_vars= {"m":"mass", "v":"volume", "rho":"density"}

and then add the properties and setters dynamically to the class (instead of explicitly):

from functools import partial

def var_getter(varname, self):
    return self._params[varname]

def var_setter(varname, self, value):
    self._params[varname] = value
    self.balance(varname)

for k,v in my_vars.items():
    setattr(Stuff, k, property(fget=partial(var_getter, v), fset=partial(var_setter, v)))

This way, you only need to write the getter and setter once.

If you want to have a couple of different getters, you can still use this technique. Store either the getter for each variable or the variables for each getter - whichever is more convenient.


Another thing that may be useful once equations get complex, is that you can keep equations as strings in your source:

density_eq = sympy.sympify(   "rho - m / v"   )

Using these two "tricks", you may even want to keep your variables and your equations defined in external text files, or maybe CSV.

Community
  • 1
  • 1
loopbackbee
  • 21,962
  • 10
  • 62
  • 97
  • Thanks, that definitely simplifies things - I'm adding some different code to my question... – flutefreak7 Mar 07 '14 at 17:15
  • Keeping things as strings makes them harder, not easier to deal with. You always have to convert them back to SymPy expressions if you want to do any kind of manipulation on them. – asmeurer Mar 08 '14 at 01:04
  • yeah I plan to use strings only when creating a new SymPy symbol, then create the expressions by combining symbols and other expressions. I plan to store the symbols and corresponding parameter strings in a dictionary together so that it's easy to either get a symbol from a parameter name for use in a solve function, or substitute a stored value to a symbol when querying the model. – flutefreak7 Mar 10 '14 at 15:48
0

Viewing your problem as a constraint-solving problem, you may also want to look at python-constraint (http://labix.org/python-constraint) in addition to SymPy.

0

Just realized that the python-constraint package doesn't apply here since the domain needs to be finite. If the domain were finite though, here's an illustrative example:

import constraint as cst

p = cst.Problem()
p.addVariables(['rho','m', 'v'], range(100))
p.addConstraint(lambda rho,m,v: rho * v == m, ['rho', 'm', 'v'])

p.addConstraint(lambda rho: rho == 2, ['rho'])     # setting inputs; for illustration only
p.addConstraint(lambda m: m == 10, ['m'])        
print p.getSolutions()

[{'m': 10, 'rho': 2, 'v': 5}]

However, since the real domain is needed here, the package doesn't apply.

bc19
  • 61
  • 5