Create Your Own Modeling Elements

The MotionSolve Python Interface allows you to create and use higher level modeling elements. The example below explains in detail. Note, the example shown below could have been done without using the MotionSolve Python Interface, however, it would have been harder to create and more difficult to use.

MotionSolve provides a rich library of general purpose and phenomena specific modeling elements. You can combine these to define new modeling elements that exhibit new behavior. This process is called aggregation. MotionSolve provides a special Composite class to represent such objects.

The Composite Class

Composite is a special class that allows you to aggregate entities, such as parts, markers, forces, differential equations, algebraic equations and other MotionSolve primitives. A composite class object behaves just like any other solver object except it is a collection of entities conveniently packaged together. They behave as atomic entities and can be instantiated multiple times. The composite class you create is inherited from the built-in composite class. The examples below shows how you can create your own composite class.

For example, assume you want to create a composite class called a low pass ButterworthFilter. The link explains the Butterworth filter is in great detail - it is summarized here. A low pass filter allows low frequency signals to pass through. It attenuates (decays) high frequency signals. A cutoff frequency defines the frequency at which the attenuation begins. As the frequency increases, the degree of attenuation increases. Butterworth filters have an order; the higher the order the more rapid is the cutoff.

The ButterworthFilter class inherits from the Composite class as shown below.

from msolve import * 

class ButterworthFilter (Composite):

Inside this class you can assign properties that you need. In this instance we wish to create a filter having an order in the range 1 ≤ order ≤ 4, with a default order of 3. Of course, you also need an input signal to operate on and an output signal to return the results. So the class with properties is:

from msolve import * 

class ButterworthFilter (Composite):
inputSignal = Function()
outputSignal = Function ()
order = Int (3)

You have to write two special methods: createChildren and updateChildren to define the behavior of this class. You can also optionally write a validate method, validateSelf.

createChildren: Used to create all the MotionSolve entities that are needed by the composite class. When the composite class object is created, the createChildren method is invoked. This happens only once, at the creation of the composite object. In createChildren you will instantiate various MotionSolve objects and define their immutable properties.
updateChildren: Used to define all the mutable properties of the class. For instance, the coefficients of the Butterworth filter depend on the order of the filter. Similarly, the input signal to act on can be changed. The updateChildren method defines the order of the filter and the input signal. They can also be defined in the createChildren method, but this is optional.
validateSelf: Checks the input to ensure that the data provided to the ButterworthFilter object is physically meaningful.

For the Butterworth filter example, the three methods may look as follows:

from msolve import * 

class ButterworthFilter (Composite):
 """The Butterworth filter is a type of signal processing filter designed to have 
 as flat a frequency response as possible in the passband. It is also referred 
 to as a maximally flat magnitude filter and it is generally used to eliminate 
 high frequency content from an input signal.
 """
inputSignal = Function ()
outputSignal = Function ()
order = Int (3)

 #- topology ----------------------------------------------------------------

def createChildren (self):
 """This method is called once when the object is created. It is used to
 define all the 'immutable' properties of the class. The Butterworth
 filter consists of the following:
 * A string containing the function expression defining the input signal
 * A VARIABLE that defines the input signal
 * A TFSISO that captures the input to output relationship
 * An ARRAY of type 'X' to hold the internal states of the TFSISO
 * An ARRAY of type 'Y' to hold the output from TFSISO
 * An ARRAY of type 'U' to hold the input signal to the TFSISO
 * A string that defines the function expression to access the output
 """
self.var = Variable (function=self.inputSignal)
x = Array (type="X")
y = Array (type="Y")
u = Array (type="U", size=1, variables=[self.var])
num = [1]
den = [1, 2, 2, 1]
self.tf = Tfsiso (label="Butterworth Filter", x=x, y=y, u=u, numerator=num,
denominator=den)
self.outputSignal = "ARYVAL({id},1)".format(id=y.id)
#- Update Children -----------------------------------------------------------
def updateChildren (self):
 """This is called when the property value changes to propagate the change
 to the child objects
 """
self.var.function = self.inputSignal
if self.order == 1:
self.tf.denominator = [1, 1]
elif self.order == 2:
self.tf.denominator = [1, math.sqrt(2), 1]
elif self.order == 3:
self.tf.denominator = [1, 2, 2, 1]
elif self.order == 4:
self.tf.denominator = [1, 2.6131, 3.4142, 2.6131, 1]
 #- validate ----------------------------------------------------------------
def validateSelf (self, validator):
validator.checkGe ("order", 1)
validator.checkLe ("order", 4)

Having created the ButterWorth composite class, you can now use it as an "atomic" object in any model. The example below is a nonlinear mass-spring-damper that oscillates in the global y-direction. Gravity is in the negative y-direction. An impulsive force of 100N is provided in the y-direction at T=1 second. The impulse ends at 1.1 seconds.

Our goal is to pass the velocity of the CM of the mass to the Butterworth filter and obtain a filtered signal as output. The response of the Butterworth filter at order=1 and order=2 are compared.

###############################################################################
## The Model ################################################################
###############################################################################
def butterworthFilter ():

m = Model ()
Units (mass="KILOGRAM",length="MILLIMETER", time="SECOND", force="NEWTON")
Accgrav (jgrav=-9800)
Integrator (error=1e-7)
ground = Part (ground=True)
                            
 # Mass
block = Part (mass=1.4702,ip=[44144.717,44144.717,73.5132,0,0,0])
block.cm = Marker (part=block, qp=[0,10,0], zp=[0,11,0])
 # Translational joint
jjm = Marker (part=ground, qp=[0,0,0], zp=[0,1,0])
jnt = Joint (type="TRANSLATIONAL", i=block.cm, j=jjm)
 # Nonlinear Spring-damper
fexp = "-10*(DY({i},{j})-10)**3  -1e-3*VR({i},{j})".format(i=block.cm.id, j=jjm.id)
nlspdp = Sforce (type="TRANSLATION", i=block.cm, j=jjm, function=fexp)

 # Impulse force
iforce = "STEP(TIME, 1, 0, 1.04, 10) *STEP(TIME, 1.05, 10, 1.1, 0)"
sfim = Marker (part=block, qp=[0,10,0], zp=[0,11,0])
impulse = Sforce (type="TRANSLATION", actiononly=True, i=sfim, j=jjm, function=iforce)

 # Filter
inputSignal = "VY({i}, {j})".format(i=block.cm.id, j=jjm.id)
m.filt = ButterworthFilter (inputSignal=inputSignal, order=1)

 # Requests
m.r1 = Request (type="DISPLACEMENT", i=block.cm, j=jjm)
m.r2 = Request (type="VELOCITY", i=block.cm, j=jjm)
m.r3 = Request (type="ACCELERATION", i=block.cm, j=jjm)
m.r4 = Request (type="FORCE", i=block.cm, j=jjm)
m.r5 = Request (f2=m.filt.inputSignal, f3=m.filt.outputSignal)
m.r6 = Request (type="FORCE", i=sfim, j=jjm)
return m
###############################################################################
## Entry Point ################################################################
###############################################################################
if __name__ == "__main__":
m = butterworthFilter ("butterworth-1")
m.simulate (type="DYNAMICS", end=5, dtout=0.005)
m.reset()
m = butterworthFilter ("butterworth-2")
m.filt.order = 2
m.simulate (type="DYNAMICS", end=5, dtout=.005)
m.reset()

Due to the action of gravity, the mass sinks down. Since there is a little damping in the system, the systems starts oscillating in the global y-direction with a decaying magnitude. At Time=1 second, the impulse force provides a jolt to the system. The system responds to the jolt and starts settling down again. The model behavior is shown below.

Figure 1.

The plot on the left below shows the velocity of the block as a function of time. This is the input to the transfer function. The plot on the right below shows the response of the order-1 and order-2 Butterworth filters. Notice how the second order filter eliminates the high frequency content.

Figure 2.

DocString

Note that Python has a powerful build in help mechanism that allows to associate documentation with modules, classes and functions. Any of such objects can be documented by adding a string constant as the first statement in the object definition. The docstring should describe what the function or class does. This allows the users and the programmer to inspect objects classes and methods at runtime by invoking the help function on a class, function or object instance.

For example the Butterworth class defined in the example above can be inspected by typing:

help (ButterworthFilter)

Provided that the module has been successfully imported in the current Python namespace, the above command will cause the Python interpreter to output the content of the docstring associated to the ButterworthFilter class. It is generally recommended to add docstrings for all objects, methods and modules in Python. The following snippet shows how module, class, method and function docstrings are defined.

"""
Assuming this is file mymodule.py, then this string, being the
first statement in the file, will become the "mymodule" module's
docstring when the file is imported.
"""
class MyClass(object):
   """The class's docstring"""
   def my_method(self):
       """The method's docstring"""
def my_function():
   """The function's docstring"""

Friction Modeling Example

Friction is a natural phenomenon that occurs in many engineering systems. Friction models contain a variety of nonlinear features such as discontinuities, hysteresis, internal dynamics, and Coulomb friction. These properties cause the friction models to be numerically stiff and computationally cumbersome. The LuGre (Lundt-Grenoble) model for friction is well known in the literature. It is capable of representing several different effects:

Dynamic friction
Static or Coulomb friction
Rate dependent friction phenomena, such as varying breakaway force and frictional lag

The model is defined in terms of the following input parameters:

VS = Stiction to dynamic friction transition velocity
MUS = Coefficient of static friction
MUD = Coefficient of dynamic friction
K0 = Bristle stiffness
K1 = Bristle damping
K2 = viscous coefficient

The friction model is formulated as follows:

V = VZ (I, J, J, J)	## Velocity in joint
P = -(v/VS)**2	## Stribeck factor
N = sqrt (FX(I, J, J)2 + FY(I,J,J)2)	## The normal force
Fc = Mud*N	## Coulomb friction
Fs = Mus*N	## Static Friction
G = (Fc + (Fs-Fc)*eP)/K0
$\dot{z}$ = v - \|v\|*z / G	## ODE defining bristle deflection z
F = -(K0z + K1zdot + K2*v)	## The friction force

The LuGre friction force consists of 2 atomic elements:

A DIFF defining the bristle deflection
A FORCE defining the friction element

The composite class can be defined as shown below.

from msolve import * 
################################################################################
## LuGre Friction Element######################################################
################################################################################
class LuGre (Composite):
"""
 Create a friction force on a translational joint using the LuGre friction model.
                            
 The LuGre friction force consists of 4 atomic elements:
 1. A DIFF defining the bristle deflection
 2. A MARKER defining the point of action of the friction force
 3. A FORCE defining the friction element
 4. A REQUEST capturing the friction force on the block
 """
joint = Reference (Joint)
vs = Double (1.E-3)
mus = Double (0.3)
mud = Double (0.2)
k0 = Double (1e5)
k1 = Double (math.sqrt(1e5))
k2 = Double (0.4)
 #- topology----------------------------------------------------------------
def createChildren (self):
 """This is called when the object is created so the children objects
 """
# The DIFF defining bristle deflection
self.diff = Diff (routine=self.lugre_diff, ic=[0,0])
# The MARKER on which the friction force acts
self.im = Marker ()
# The FORCE defining the friction force
self.friction = Sforce (type="TRANSLATION", actiononly=True,
                                routine=self.lugre_force)
# The REQUEST capturing the friction force
self.request = Request (type="FORCE", comment="Friction force")
 #- Update Children -----------------------------------------------------------
def updateChildren (self):
 """This is called when the property value changes to propagate the change
 to the child objects
"""
self.friction.i = self.joint.i
self.friction.j = self.joint.j
self.im.setValues (body=self.joint.i.body, qp=self.joint.i.qp, zp=self.joint.i.zp)
self.request.setValues (i=self.im, j=self.joint.j, rm=self.joint.j)
 #- validate ----------------------------------------------------------------
def validateSelf (self, validator):
validator.checkGe0 ("VS")
validator.checkGe0 ("MUS")
validator.checkGe0 ("MUD")
validator.checkGe0 ("K0")
validator.checkGe0 ("K1")
validator.checkGe0 ("K2")
if self.mud > self.mus:
msg = tr("Mu dynamic({0}) must be <= Mu static ({1})",
         self.mud, self.mus)
         validator.sendError(msg)
 #- Lugre Diff ------------------------------------------------------------
def lugre_diff (self, id, time, par, npar, dflag, iflag):
 "Diff user function"
i = self.joint.i
j = self.joint.j
vs = self.vs
mus = self.mus
mud = self.mud
k0 = self.k0
z = DIF(self)
v = VZ(i,j,j,j)
N = math.sqrt (FX(i,j,j)**2 + FY(i,j,j)**2)
fs = mus*N
fc = mud*N
p = -(v/vs)**2
g = (fc + (fs - fc) * math .exp(p))/k0
if iflag or math.fabs(g) <1e-8:
return v
else:
return v - math.fabs(v) * z / g
 #- friction_force -----------------------------------------------------------
def lugre_force (self, id, time, par, npar, dflag, iflag):
 "Friction Force user function"
i = self.joint.i
j = self.joint.j
diff = self.diff
k0 = self.k0
k1 = self.k1
k2 = self.k2
v = VZ (i,j,j,j)
z = DIF (diff)
zdot = DIF1 (diff)
F = k0*z + k1*zdot + k2*v
return -F