Programmatic construction of the three-link pendulum

From neuromechanicwiki
Jump to: navigation, search

In the previous tutorial we added a neural controller in Python to the muscles of a three-link pendulum. By this point you may have noticed that it is pretty tedious to build a model by writing the XML body file by hand. It's especially tedious

  1. If you have a form that you basically want to build up from a few simple elements. In our case we only have a few muscle forms and a couple body forms that we are essentially transforming many times.
  2. If you already have a Neuromechanic body file and you just want to tweak a few elements many times (e.g. parameters search)
  3. If you have a model in some other defined structure (e.g. opensim file) and you just want to translate it to the Neuromechanic body file.

For this reason it is useful to have a programmatic way to construct a body file. In this tutorial we'll go over how to do that in the Neuromechanic engine with Python.

You will need...

We need some tools to assist in constructing the file.

Download the following files:

  • nmextension.py uses Python dictionaries and the xml.ElementTree module to convert back and forth between an xml file and a special Python object based on Python dictionaries.
  • pendulum3control.py this is the controller module that we build in the last tutorial


Save these files into a common directory that is easily accessible.

A simple example to begin with

The following code will create a single bodied Neuromechanic file, load it in the engine and display it. Copy the text into a *.py module and import the module from the command line.

# First import the modules we'll need
import os
import numpy as np
import neuromechanic as nm
import nmextension as nme
 
x = nme.Build()
 
# Add the root body at the global origin
x.addRigidBody('root','ground',FrameLocation=[0,0,0],CenterOfMass=[0,0.12,0],InertialEllipsoid=[.04,.2,.04],Density=1000)
# The root body will have 6 degrees of freedom. All locked
x.addDoF(-1,'tx','Translation',Axis=[1.,0.,0.],Locked='true')
x.addDoF(-1,'ty','Translation',Axis=[0.,1.,0.],Locked='true')
x.addDoF(-1,'tz','Translation',Axis=[0.,0.,1.],Locked='true')
x.addDoF(-1,'rx','Rotation',Axis=[1.,0.,0.],Locked='true')
x.addDoF(-1,'ry','Rotation',Axis=[0.,1.,0.],Locked='true')
x.addDoF(-1,'rz','Rotation',Axis=[0.,0.,1.],Locked='false')
 
# Get the current working directory, this is where the file will be saved
pth = os.getcwd();
 
# Create the xml file    
err = x.write(os.path.join(pth,'testpendulum.nmcb'))
print('If the error ('+str(err)+'==0) then the file was successfully written to disk')
 
# Load the file into the engine
nm.file_load(os.path.join(pth,'testpendulum.nmcb'))
 
# Draw it
nm.gui_renderall()

Constructing a model body file in Neuromechanic with Python

The following instructions will go step-by-step in building a python file that will programmatically build a Neuromechanic model file. You can either:

or

  • Build the module from scratch by creating a new text file named constructpendulum.py.


The first thing to do is to import the modules we will need. Add these lines to the top of the file:

  1. import os
  2. import numpy as np
  3. import neuromechanic as nm
  4. import nmextension as nme

Each import statement loads a set of library functions that will be associated with a shortened name specified by the as keyword. For example, we will be accessing neuromechanic using the shortened name nm.

Next we define a function that will add the additional element our Monoarticular motor neurons need

  1. def mononeuron(x,ind,mname,phi):
  2.     '''
  3.     This will add some elements to the "ind"th Monoarticular Motor Neuron for use with
  4.     the python controller. The only inputs beside the structure "x" and the index
  5.     of the neuron to be added to "ind" are the name of the muscle to reference
  6.     "mname" and the phase "phi"
  7.     <Neuron Name="MNLink1B" Type="NeuralNet">
  8.       ...
  9.       <drive>0.</drive>
  10.       <phase>phi</phase>
  11.       <musclenode DOMNode=".^.^.^.Muscles.Muscle{Name=mname}"/>
  12.       <FunctionCall Type="Simulation" Name="MonoarticularActivation">
  13.         <Arguments>musclenode phase drive</Arguments>
  14.         <Return>NeuronOutput</Return>
  15.       </FunctionCall>
  16.       <FunctionCall Type="Equilibration" Name="MonoarticularActivation">
  17.         <Arguments>musclenode phase drive NeuronOutput</Arguments>
  18.         <Return>NeuronOutput drive</Return>
  19.       </FunctionCall>
  20.       <NeuronOutput>1.</NeuronOutput>
  21.     </Neuron>
  22.     '''
  23.  
  24.     x.NeuromechanicFile.Neurons.Neuron[ind]['drive'].x_ = 0.
  25.     x.NeuromechanicFile.Neurons.Neuron[ind]['phase'].x_ = phi
  26.     x.NeuromechanicFile.Neurons.Neuron[ind]['musclenode'] = xd.XmlDictObject()
  27.     x.NeuromechanicFile.Neurons.Neuron[ind].musclenode['DOMNode_'] = '.^.^.^.Muscles.Muscle{Name='+mname+'}'
  28.     x.NeuromechanicFile.Neurons.Neuron[ind]['FunctionCall'] = []
  29.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall.append(xd.XmlDictObject())
  30.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[0]['Type_'] = 'Simulation'
  31.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[0]['Name_'] = 'MonoarticularActivation'
  32.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[0]['Arguments'] = 'musclenode phase drive'
  33.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[0]['Return'] = 'NeuronOutput'
  34.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall.append(xd.XmlDictObject())
  35.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[1]['Type_'] = 'Equilibration'
  36.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[1]['Name_'] = 'MonoarticularActivation'
  37.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[1]['Arguments'] = 'musclenode phase drive NeuronOutput'
  38.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[1]['Return'] = 'NeuronOutput drive'
  39.  
  40.     return x

Next we define a function that will add the additional element our Biarticular motor neurons need

  1. def bineuron(x,ind,mname):
  2.     '''
  3.     This will add some elements to the "ind"th Biarticular Motor Neuron for use with
  4.     the python controller. The only inputs beside the structure "x" and the index
  5.     of the neuron to be added to "ind" are the name of the muscle to reference
  6.     "mname" and the phase "phi"
  7.     <Neuron Name="MNLink1B" Type="NeuralNet">
  8.       ...
  9.       <drive>0.</drive>
  10.       <phase>phi</phase>
  11.       <musclenode DOMNode=".^.^.^.Muscles.Muscle{Name=mname}"/>
  12.       <FunctionCall Type="Simulation" Name="MonoarticularActivation">
  13.         <Arguments>musclenode phase drive</Arguments>
  14.         <Return>NeuronOutput</Return>
  15.       </FunctionCall>
  16.       <FunctionCall Type="Equilibration" Name="MonoarticularActivation">
  17.         <Arguments>musclenode phase drive NeuronOutput</Arguments>
  18.         <Return>NeuronOutput drive</Return>
  19.       </FunctionCall>
  20.       <NeuronOutput>1.</NeuronOutput>
  21.     </Neuron>
  22.     '''
  23.  
  24.     x.NeuromechanicFile.Neurons.Neuron[ind]['drive'] = 0.
  25.     x.NeuromechanicFile.Neurons.Neuron[ind]['musclenode'] = xd.XmlDictObject()
  26.     x.NeuromechanicFile.Neurons.Neuron[ind].musclenode['DOMNode_'] = '.^.^.^.Muscles.Muscle{Name='+mname+'}'
  27.     x.NeuromechanicFile.Neurons.Neuron[ind]['FunctionCall'] = []
  28.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall.append(xd.XmlDictObject())
  29.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[0]['Type_'] = 'Simulation'
  30.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[0]['Name_'] = 'BiarticularActivation'
  31.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[0]['Arguments'] = 'musclenode drive'
  32.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[0]['Return'] = 'NeuronOutput'
  33.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall.append(xd.XmlDictObject())
  34.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[1]['Type_'] = 'Equilibration'
  35.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[1]['Name_'] = 'BiarticularActivation'
  36.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[1]['Arguments'] = 'musclenode drive NeuronOutput'
  37.     x.NeuromechanicFile.Neurons.Neuron[ind].FunctionCall[1]['Return'] = 'NeuronOutput drive'
  38.  
  39.     return x

Next we'll define the key dimensions and other parameters of our three-link pendulum model. We'll use those parameters to help define a few other useful variables.

  1. # Parameters *********************************************************
  2. fname = 'junk.nmcb'                 # The name of the file to save
  3. hingeaxis = np.array([0.,0.,1.])    # Direction of hinge axis
  4. longaxis = np.array([0.,1.,0.])     # Direction of long axis of pendulum (perpendicular to hingeaxis)
  5. jointdistance = .24;                # Distance between joints
  6. segmentlength = .2;                 # Major axis of equivalent inertial ellipsoid
  7. segmentdiameter = .04;              # Length of the other axes of equivalent inertial ellipsoids
  8. musclespace = .03;                  # Attachment distance of muscles from joint
  9. maxforce = 80                       # Maximum isometric muscle force
  10. numberoflinks = 3                   # Number of links to add
  11. # Parameters *********************************************************
  12.  
  13. # The third axis is perpendicular to the hinge and long axis
  14. third = np.abs(np.cross(longaxis,hingeaxis))
  15.  
  16. # Make two arrays to hold equivalent inertial ellipsoid dimensions
  17. # The first is for the root body and should be a sphere with a diameter of "segmentdiameter"
  18. eieroot = (third+hingeaxis+longaxis)*segmentdiameter
  19. # The second is for the pendulums, the ellipsoid is "segmentlength" long in the "longaxis"
  20. # direction and "segmentdiameter" wide in the other two dimensions
  21. eie = (third+hingeaxis)*segmentdiameter + longaxis*segmentlength;
  22.  
  23. # Make a muscle attachment point 
  24. musattach = musclespace*third

The next step is to initialize a Python object that will hold all the data to be written to the xml file. Here is where we will start using nmbuild.py so now is a good time to review some of those functions. Essentially, nmbuild.py contains helper routines to populate the structure.

  1. # Initialize the structure to hold the "xmldict" object
  2. x = nme.Build()
  3.  
  4. # Add the root body at the global origin
  5. x.addRigidBody('root','ground',[0,0,0],[0,0,0],eieroot)
  6. # The root body will have 6 degrees of freedom. All locked
  7. x.addDoF(-1,'tx','Translation',[1.,0.,0.],Locked='true')
  8. x.addDoF(-1,'ty','Translation',[0.,1.,0.],Locked='true')
  9. x.addDoF(-1,'tz','Translation',[0.,0.,1.],Locked='true')
  10. x.addDoF(-1,'rx','Rotation',[1.,0.,0.],Locked='true')
  11. x.addDoF(-1,'ry','Rotation',[0.,1.,0.],Locked='true')
  12. x.addDoF(-1,'rz','Rotation',[0.,0.,1.],Locked='true')

Next we're going to start the loop that adds the links of the pendulum.

  1. # Initialize some variables to assist in building the links
  2. grandparent = [];
  3. grandfloc = [];
  4. phase = 0
  5. # Loop through the number of links, build a link each time
  6. for jj in range(numberoflinks):
  7.  
  8.     # Increment the name of the current body we're building
  9.     name = 'link' + str(jj+1)
  10.  
  11.     # The FrameLocation body and position depend on whether this is
  12.     # the first or a subsequent link
  13.     if jj==0:
  14.         # The FrameLocation of the first body is at the orgin of "root"
  15.         floc = np.array([0,0,0])
  16.         parent = 'root'
  17.     else:
  18.         # Otherwise we use our parameters to get the FrameLocation
  19.         floc = jointdistance*np.array(longaxis)
  20.         parent = 'link' + str(jj)
  21.  
  22.     # The center of mass is
  23.     com = longaxis*jointdistance/2.
  24.  
  25.     # Add a new body
  26.     x = nmb.addrigidbody(x,name,parent,floc,com,eie)
  27.  
  28.     # Create a rotational <DegreeOfFreedom> for this body in the "hingeaxis" direction
  29.     x = nmb.adddof(x,-1,name+'rz','Rotation',hingeaxis,State=[0,0])

At each iteration of the loop we will also add monoarticular muscles and motor neurons to the new body.

  1.     # Add a monoarticular pair of muscles
  2.     # Lets make tendon slack length equal to the attachment offset, just because
  3.     tsl = musclespace           
  4.     # Guess what the muscle length will be and subtract the tendon length to get optimal fiber length
  5.     lf0 = jointdistance/2.-tsl  
  6.  
  7.     # The first muscle will originate on the right and insert at the center of mass
  8.     path = [{'Parent':parent,'x':floc+musattach},{'Parent':name,'x':com}]
  9.     x = nmb.addzajacmuscle(x,name+'R',path,'MN'+name+'R',maxforce,lf0,tsl)
  10.     x = nmb.addneuron(x,'MN'+name+'R',EqCCost=-.5)
  11.     # The Second muscle will originate on the left and insert at the center of mass
  12.     path = [{'Parent':parent,'x':floc-musattach},{'Parent':name,'x':com}]
  13.     x = nmb.addzajacmuscle(x,name+'L',path,'MN'+name+'L',maxforce,lf0,tsl)
  14.     x = nmb.addneuron(x,'MN'+name+'L',EqCCost=-.5)
  15.  
  16.     # These functions add the extra variables to the monoarticular muscles
  17.     # the extra variables are used in the python controller
  18.     # One of the variables is phase which we will increment by 180 degrees
  19.     # at each link and the left and right muscles are 180 degrees out of phase
  20.     phase = phase + 180
  21.     x = mononeuron(x,-2,name + 'R',phase)
  22.     x = mononeuron(x,-1,name + 'L',phase+180)

The last thing that is done each iteration of the loop is biarticular muscles are added (as long as this is not the first link.

  1.     # Add a biarticular pair of muscles if this is not the first link
  2.     if grandfloc != []:
  3.         # Lets make tendon slack length equal to the attachment offset, just because
  4.         tsl = musclespace
  5.         # Guess what the muscle length will be and subtract the tendon length to get optimal fiber length
  6.         lf0 = jointdistance-tsl
  7.  
  8.         # The first muscle will originate on the right of the "grandparent" body 
  9.         # and insert on the right of the current body
  10.         path = [{'Parent':grandparent,'x':grandfloc+musattach},{'Parent':name,'x':np.array([0.,0.,0.])+musattach}]
  11.         x = nmb.addzajacmuscle(x,name+'BR',path,'MN'+name+'BR',maxforce,lf0,tsl)
  12.         x = nmb.addneuron(x,'MN'+name+'BR')
  13.         # The second muscle will originate on the left of the "grandparent" body 
  14.         # and insert on the left of the current body
  15.         path = [{'Parent':grandparent,'x':grandfloc-musattach},{'Parent':name,'x':np.array([0.,0.,0.])-musattach}]
  16.         x = nmb.addzajacmuscle(x,name+'BL',path,'MN'+name+'BL',maxforce,lf0,tsl)
  17.         x = nmb.addneuron(x,'MN'+name+'BL')
  18.  
  19.         # These functions add the extra variables to the monoarticular muscles
  20.         # the extra variables are used in the python controller
  21.         x = bineuron(x,-2,name + 'BR')
  22.         x = bineuron(x,-1,name + 'BL')
  23.  
  24.     # Store the name of the parent for this body (will be grandparent for the next body)
  25.     grandfloc = floc;
  26.     grandparent = parent;

After all the links and muscles are added we need to add some additional details about the model

  1. # Add in "afl" (active force length), "pfl" (passive force length), and "afv" 
  2. # (active force velocity) splines.
  3. x = nmb.addmusclesplines(x);
  4.  
  5. # Change equilibration cost function to MinimumForce
  6. x.NeuromechanicFile.Parameters['EquilibrationType'] = 'MinimumForce'
  7.  
  8. # Add gravity vector
  9. x = nmb.addgravity(x,[0.,-9.80665,0.])
  10.  
  11. # Add python modules and routines for controller
  12. x = nmb.addpythonmodule(x,[{'Module':'pendulum3control','Routine':['BiarticularActivation','MonoarticularActivation']}])

Finally, we can build, load, and render the new file.

  1. # Get the current working directory, this is where the file will be saved
  2. pth = os.getcwd();
  3.  
  4. # Create the xml file    
  5. xd.ConvertDictToXml(x,os.path.join(pth,fname))
  6.  
  7. # Load the file into the engine
  8. nm.system_load(os.path.join(pth,fname))
  9.  
  10. # Draw it
  11. nm.gui_renderall()

So how should we run this? Make sure this constructpendulum.py is saved and is in the same directory as pendulum3control.py. Open Neuromechanic and type each of the following lines into the command line. NOTE in windows you will need to escape backslashes in the directory (e.g. D:\\Users\\...)

  1. pth = '...path to directory containing this constructpendulum.py module';
  2. import sys
  3. import os
  4. os.chdir(pth)
  5. sys.path.append(pth)
  6. import constructpendulum

Did it work? Unfortunately right now the gui doesn't understand when a model has been loaded from python. To actually run a simulation you will need to close the model you just made and reopen it from the gui. Try running a simulation to make sure the model behaves as it did in the last tutorial. Now, tweak some of the parameters (e.g. segmentlength) and see what happens!

In the next tutorial, Extending the model to a quadruped, all we're going to do is add a rotational degree of freedom to each link, twice as many muscles per pendulum and wrap that all in one more loop that will replicate and transform the three-links four times to make a quadruped.