# Neural control of a three-link pendulum

{{#ev:youtube|lS4zox1y904|400|right|Wiggling Pendulum}} In the last tutorial we added mono and biarticular muscles to a model of a three-link pendulum. We used a constant "excitation" level to the muscles and learned how to set that constant value to balance the external forces acting on the links of the pendulum, thereby equilibrating the model. In this tutorial, we will implement some simple control principles into the model. "Neural control" is used somewhat loosely here: we will implement simple linear control gains on muscle length and velocity representing neural control based on muscle spindle sensory feedback.

In this tutorial we will add a phasic activation pattern to the monoarticular muscles to make the inverted pendulum wiggle. as shown in the movie on the right. By the end of this tutorial you will have made two files which, if you want to cheat, can just be downloaded here: pendulum3control.nmcb and pendulum3control.py

A second implementation of this controller can be found in pendulum3control_faster.nmcb. This is a much faster implementation since it only calls python from the engine once per integration step. We don't spend any time in the tutorial going over the differences between the two implementations but hopefully they will be apparent by comparing the two files.

Both nmcb files produce the wiggling pendulum shown in the movie

## Overview of scripted control in Neuromechanic

We are finally getting to one of the primary motivations for developing Neuromechanic: to have a simple but powerful way to try out different control strategies with biomechanical models. Up until this point we haven't needed to know much about programming, an understanding of XML syntax has been sufficient. Some very simple and low-level control principles as well as neural networks can be implemented through the <Synapse> elements in the XML body file. However, for more sophisticated control, we need:

1. Access to much more information about the model than can be provided in <Synapse> elements. You can think of this as a richer sensory data set if you like.
2. The ability to use that information in any way we please in order to set the excitation of the muscles.

In Neuromechanic, many of the parameters and variables which describe the anatomy and dynamics of an animal model are stored in a model tree which reflects the structure of the XML body file. Access to the model tree is granted through a C++ and Python application programming interface (API). The API also provides access to many of the core Neuromechanic routines for running a simulation. The Neuromechanic GUI uses the C++ API to interface with the Neuromechanic engine.

This tutorial makes use of the Python version of the API. Python has a relatively low barrier to entry if you're new to programming or have only been programming with Matlab. I would like to believe that you will be able to comprehend this tutorial and the next few even if you don't know Python. I'm going to try to comment my code in a way that I think I would've liked to see the first Python routines I looked at.

Some excellent tutorials exist on the web for python and a quick search should present you with an abundance of sources. A good beginners guide is offered through some Google python tutorials.

## Two ways to use Python in Neuromechanic

The embedded Python interpreter can be used in two ways.

1. In this tutorial we will write Python helper functions for evaluating <Neuron> output. The Neuromechanic engine is the master and is responsible for calling these functions when it is appropriate.
2. In the next tutorial Python will be the master. We will write a Python function that calls various routines in the Neuromechanic engine. It is initiated by the user from the Neuromechanic command line.

## The Neuromechanic Model Tree

Before we get into the details of the Python functions further explanation about the Neuromechanic model tree is needed. The model tree is similar to a document object model DOM. The Neuromechanic model tree isn't a real DOM, that's just the closest thing to it. The model tree is really just the representation of the (extended) XML file in memory. Here are the key functional differences between a traditional DOM and Neuromechanic's model tree:

### Text Nodes are not text nodes

First, the "Text Node" in a traditional DOM is replaced in the Neuromechanic model tree with either text, a real matrix, or an integer matrix, depending on the content of the node. In this code for example,

1. <Random name="Bob" Age="4">
2.    <Thing>
3.       1 0 0
4.       0 1 0
5.       0 0 1
6.    </Thing>
7. </Random>

the Neuromechanic XML parser would recognize that the content of <Thing> is a 3x3 integer matrix. If you queried the value of <Thing> in the model tree you would not find a character string, but a 3x3 integer matrix. If <Thing> was something that Neuromechanic recognized as a particular part of a body file that needed to be a real matrix then it would automatically convert <Thing> from integer to real.

An exception to this number recognition is that attributes can ONLY be character strings. For the example shown, the value of the attribute "Age" would be represented in the model tree with a character string "4" rather than the number 4.

The bottom line is: when you query the model tree about a Node, you don't just get text

### Reserved Elements

Neuromechanic expects certain content in an XML body file. Usually, that content must be supplied by the user. However, there are some elements that can safely be generated by Neuromechanic automatically and placed in the model tree when the XML file is parsed. For example, if the <View> element is missing from a body file Neuromechanic will automatically generate all that content. When the model tree is printed to an XML file, that auto-generated content will appear in the new XML file.

Some elements that are auto-generated in the model tree SHOULD NOT be supplied by the user and will usually not appear in the XML generated from the model tree. These are the Reserved elements. An example of a reserved element is muscle force. We have talked about <Muscle> elements before and you will recall that nowhere in the element do we specify a child called <MuscleForce> to hold the force generated in the muscle. It's fairly obvious why this is not done: Muscle force is not an independent parameter. It is a variable that is dependent on many independent parameters, model states, time, and the particular muscle model in use. In other words it is not specified directly by the user. However, muscle force is obviously pretty important and something a user would want access to. In fact <MuscleForce> IS a reserved element and thus does exist in the model tree and can be accessed by the user, even though it should never be changed by the user. If the user does specify something in a <MuscleForce> element it will be overwritten when the file is parsed. Two additional reserved elements in each <Muscle> are the <MusculotendonLength> and <MusculotendonVelocity>.

The bottom line is: everything in the XML gets put in the model tree but the model tree does not contain ONLY the information that you can see in the XML

### Additional contents of nodes

Besides the text nodes each node has additional information content. For example, you can query each node about whether it is "printable". The reserved elements discussed previously are not printable and therefore they don't show up in the XML.

The bottom line is: sometimes we need to get information about a node besides its value

## Initializing Callable Python Modules and Functions

We'll come back to all this very soon but we're going to switch gears for a second.

A Python module is just a file that contains Python code. A Python module can contain one or more (or less) functions. Python functions are just bundles of code that go nicely together. In this case we'll make a Python function to control the activation of our monoarticular muscles (let's call it MonoarticularActivation). We'll put this function into a single file (module) named "pendulum3control.py". That module and the contained functions need to be initialized before they can be used in Neuromechanic. This is done by creating a <PythonModule> element in the <Resources> element. module name is set as the attribute "Name" in the element <PythonModule> and each function contained in the module that needs to be called from Neuromechanic are included in <PythonRoutine> elements.

Open your file from the last tutorial and add the following code to the resources section

1. <Resources>
2.    ...
3.     <PythonModule Name="pendulum3control">
4.       <PythonRoutine>MonoarticularActivation</PythonRoutine>
5.     </PythonModule>
6.    ...
7. </Resources>

WARNING: Usually Neuromechanic is case-insensitive but when dealing with python, you should be careful that case is matched exactly

## Interfacing Callable Python Functions from Neurons

Next we need to make the interface for the functions. The interface to the function is a <FunctionCall> element that is placed in the <Neuron> element which calls the function.

1.     <Neuron Name="MNLink1A" Type="NeuralNet">
2.       ...
3.       <FunctionCall Name="MonoarticularActivation"/>
4.     </Neuron>

<FunctionCall> contains the name of the function in the Name attribute (which, of course should match the name initialized in <PythonRoutine>)

1.     <Neuron Name="MNLink1A" Type="NeuralNet">
2.       <!-- These are elements that we've talked about previously, each tag is a reserved keyword.
3.        Notice that we have synapse inputs to the Neuron from the MusculotendonLength, MusculotendonVelocity,
4.        and MuscleForce of "MuscleLink1A".  However only the gain of MusculotendonLength is non-zero here so
5.        it alone of the three will affect the neural output.-->
6.       <EqCCost>-5.E-1</EqCCost>
7.       <EqBounds>1.E-2 9.5E-1</EqBounds>
8.       <OutputBounds>0. 1.</OutputBounds>
9.       <Synapse Name="whatever" Type="constant">-9.395453501482 0.</Synapse>
10.       <Synapse Name="MuscleLink1A" Type="MusculotendonLength">8.E1 0.</Synapse>
11.       <Synapse Name="MuscleLink1A" Type="MusculotendonVelocity">0. 0.</Synapse>
12.       <Synapse Name="MuscleLink1A" Type="MuscleForce">0. 0.</Synapse>
13.       <!-- NeuronOutput is a reserved element and will ALWAYS be created and set to a
14.       real scalar when the file is parsed. If the user sets this value in the XML file
15.       it will likely be overwritten when the file is parsed. However, it CAN be set
16.       by the user in the python function we are building for this neuron.-->
17.       <NeuronOutput>5.000000000004E-1</NeuronOutput>
18.       <!-- You can create new elements as long as they don't collide with other
19.       names in <Neuron> (e.g. Synapse). We recommend that to do this you use a leading
20.       underscore in the tag name. DO NOT use a trailing underscore. The user-defined elements
21.       can hold character strings, integer or real scalars or integer or real matrices, just
22.       like the rest of the nmcb file-->
23.       <_phase>0.</_phase>
24.       <!-- The python definition for this function would look something like this
25.       def MonoarticularActivation(neuronp):
26.          ...
27.          return-->
28.       <FunctionCall Name="MonoarticularActivation"/>
29.     </Neuron>

Remember, the goal is to add a phasic component to the activation of the Monoarticular muscles. The <_phase> element contains information that will be used in the phasic controller we are about to build.

## Updating the <Neuron> elements to interface with Python

Here is what <Neuron> should look like for each of the six monoarticular muscles

    <Neuron Name="MNLink1A" Type="NeuralNet">
<EqCCost>-5.E-1</EqCCost>
<EqBounds>1.E-2 9.5E-1</EqBounds>
<OutputBounds>0. 1.</OutputBounds>
<Synapse Name="whatever" Type="constant">0. 0.</Synapse>
<Synapse Name="MuscleLink1A" Type="MusculotendonLength">8.E1 0.</Synapse>
<Synapse Name="MuscleLink1A" Type="MusculotendonVelocity">0. 0.</Synapse>
<Synapse Name="MuscleLink1A" Type="MuscleForce">0. 0.</Synapse>
<_phase>0.</_phase>
<FunctionCall Name="MonoarticularActivation"/>
<NeuronOutput>0.</NeuronOutput>
</Neuron>

Copy the contents to each mono-articular muscle but make sure that

1. The name of each <Neuron> is preserved and
2. The <Synapse> element's name attributes are set to the appropriate muscle.

The biarticular muscles will be controlled slightly differently. We'll give them synaptic feedback from the homonymous muscle's length and velocity (force gain is set to zero). They WILL not be controlled by a python function and should look like this.

    <Neuron Name="MNLink12A" Type="NeuralNet">
<EqCCost>0.</EqCCost>
<EqBounds>1.E-2 9.5E-1</EqBounds>
<OutputBounds>0. 1.</OutputBounds>
<Synapse Name="whatever" Type="constant">-1.189906 0.</Synapse>
<Synapse Name="MuscleLink12A" Type="MusculotendonLength">5. 0.</Synapse>
<Synapse Name="MuscleLink12A" Type="MusculotendonVelocity">2.E-1 0.</Synapse>
<Synapse Name="MuscleLink12A" Type="MuscleForce">0. 0.</Synapse>
<NeuronOutput>0</NeuronOutput>
</Neuron>

Copy the contents to each biarticular muscle but again, make sure that

1. The name of each <Neuron> is preserved and
2. The <Synapse> element's name attributes are set to the appropriate muscle.

## Python Calls

Alright, the XML body file is set up, now all we have to do is make the MonoarticularActivation function. In the <Resources> section we said that it belongs to a module called pendulum3control. To create that module we need to create a file named pendulum3control.py copy the following to the file

1. # Our module makes use of the code in other modules. We need to "import" those modules here
2. 
3. # - math has some functions and variables (e.g. sin, pi) that will be useful in our controller
4. # when you import a module you create a namespace associated with it. That means if we want
5. # to call some code in math (e.g. sin()) then we'd call it like this  y = math.sin(x)
6. import math
7. 
8. # - nmextenstion is a module that contains the neuromechanic api and additional useful code
9. # by adding the "as nme" in the import statement we rename the namespace. So if we want
10. # to use some code in nmextension (e.g. The class "Tree") instead of calling nmextension.Tree(...)
11. # we will call nme.Tree(...). If you opened up the file nmextension.py you would see import
12. # statements at the top of it just like in this file. One of them is "import neuromechanic as nm"
13. # the module "neuromechanic" is the neuromechanic api, it's built-into the engine and a
14. # description of its contents can be found on the wiki (wiki.neuromechanic.com). So since it
15. # was imported "as nm" by nmextension and nmextension is imported here "as nme" that means you
16. # can access the neuromechanic contents (for example the function version()) with y = nme.nm.version()
17. import nmextension as nme
18. 
19. 
20. # nmtree will be an "instance" of a "class" called Tree that is stored in nmextension.
21. # We initially set it to "None". Which basically means it starts empty.
22. nmtree = None
23. 
24. def MonoarticularActivation(neuronp):
25.     '''
26.     Define a function "MonoarticularActivation" with a single input "neuronp". This single
27.     input is the format for ALL functions called by a <Neuron> in Neuromechanic. "neuronp" is
28.     a pointer to the <Neuron> element itself.
29. 
30.     For those unfamiliar with python syntax, notice the colon after the def statement. Everything
31.     that is a part of MonoarticularActivation should have indententation after that colon.
32.     There is no "end" statement for python functions, The function includes everything below
33.     the def statement that maintains this indentation.
34. 
35.     The function will return an activation value (NeuronOutput)
36.     Each function called by a <Neuron> in Neuromechanic is passed a pointer to that <Neuron>.
37.     (NeuronOutput) as well as a descending drive (drive) value.
38.     '''
39. 
40.     # nmtree is declared in the module outside this function, by declaring it global
41.     # here we can set it's value from inside this function and its value will be preserved
42.     # between successive calls to MonoarticularActivation or BiarticularActivation
43.     global nmtree
44. 
45.     if nmtree is None:
46.         # Tree is a special class in the nmextension module that lets us interface with
47.         # the parts of the model tree. If tree is called with the keyword "model" then
48.         # it will create a map to the entire model tree. "nmtree" will be our instance
49.         # of that map. This interfacing only needs to be done once so we place it behind
50.         # this conditional statement, once it is set it will no longer be "None" so this
51.         # interfacing only happens once.
52.         nmtree = nme.Tree('model')
53. 
54.     # If you pass the "Tree" class a pointer to a location in the model tree then instead
55.     # of mapping the entire tree it maps the tree starting from the location you passed it.
56.     # Here we pass Tree a pointer to the <Neuron> element that called this function.
57.     # So the all descendants of that calling <Neuron> will be mapped in the variable "neuron"
58.     neuron = nme.Tree(neuronp)
59. 
60.     # Now that we have a map to the <Neuron> that called this function we can retrieve some
61.     # other parameters stored there the <_phase> and <_drive> elements.
62.     # ... = neuron._phase will get us to a pointer to the _phase node. To actually retrieve
63.     # the contents of the node we append ".x_"
64.     phase = neuron._phase.x_            # radians.    Phase shift of sinusoid
65. 
66.     # Two additional parameters we need are <_frequency> and <_amplitude>. That were stored
67.     # in the <Neurons> element. We'll get them by navigating "nmtree"
68.     freq = nmtree.NeuromechanicFile.Neurons._frequency.x_   # Hz            The frequency of the oscillating activation signal
69.     amp = nmtree.NeuromechanicFile.Neurons._amplitude.x_    # unitless      Amplitude of the oscillating activation signal
70.     # We'll also want the current simulation time which is stored in the reserved element
71.     # <Dynamic>
72.     t = nmtree.NeuromechanicFile.Dynamic.Time.x_
73. 
74.     # Calculate the oscillating activation signal
75.     sinudrive = amp*math.sin(2.*math.pi*freq*t + phase*math.pi/180.)
76. 
77.     # Calculate the output value
78.     neuron.NeuronOutput.x_ = neuron.NeuronOutput.x_ + sinudrive
79. 
80.     # Python has several ways to group data together (like matlab has structures,
81.     # arrays and cell arrays. While there is no direct equivalency it seems to me
82.     # that cell arrays in matlab are closest to lists (e.g. [1,2,'hi']) and tuples
83.     # (e.g. (1,2,'hi')) in python. Structures in matlab are closest to dictionaries
84.     # in python. Python has no restrictions regarding what kind of data is returned
85.     # from a function but Neuromechanic expects a tuple or a single value, not a
86.     # list, if you want to return a list then put it in a tuple before returning it.
87.     return

To run a simulation with this controller make sure the python module and body file are in the same directory. Load the file and do the following:

• Initial posture is zero for all joint angles (an upright inverted pendulum).
• Make sure the <_phase> is 0 degrees for motor neurons MNLink1A, MNLink2B, and MNLink3A and 180 degrees for Neurons MNLink1B, MNLink2A, and MNLink3B.
• Set <EqCCost> of monoarticular muscles to -0.5 and biarticular muscles to 0
• Changed <EquilibrationType> to minimumForce
• Changed <EndTime> to 3 seconds
• Equilibrate
• Begin Simulation

## Speed Improvement

Ok, that's just about it. You've seen how to format your nmcb file and create a python module to make your own Neural controller. Before we finish I just want to add one more wrinkle. The implementation we just built makes a call from the Neuromechanic engine to the Python routine 6 times (once fore each monoarticular muscle) per integration step. There is a performance penalty for just calling Python, and even more so when you are repeating some calculations every time you call it. Another way we could structure our controller is to make one <Neuron> the "Controller" and have the motor neurons receive synaptic input from that controller. In this way you'd only call Python once per integration step. We won't go through the details here but you can download this more efficient implementation here: pendulum3control_faster.nmcb. Here is the controller for this nmcb file. You should copy it to the same module (pendulum3control.py) you just created.

1. def MonoarticularController(neuronp):
2.     '''
3.     Here is another function that is not called by pendulum3control.nmcb. This is for use
4.     with the function pendulum3control_faster.nmcb. One problem with what we just learned
5.     is that each motor neuron is calling the python function and repeating several calculations
6.     This can chew up a lot of simulation time. Instead in pendulum3control_faster.nmcb we make
7.     a single <Neuron> (call it "controller"). "Controller" sets the value for the sinusoidal drive
8.     as parameters within the "Controller" <Neuron>. For example the sinusoidal component of the drive
9.     to motor neuron "_MNLink1A" is set to the internal component <_MNLink1A_act> inside the controller
10.     neuron. Then another synapse is added to the motor neuron "_MNLink1A" that adds this value to the
11.     motor neuron output without calling python again.
12.     '''
13.     global nmtree
14.     neuron = nme.Tree(neuronp)
15.     if nmtree is None:
16.         nmtree = nme.Tree('model')
17. 
18.     t = nmtree.NeuromechanicFile.Dynamic.Time.x_ # s    Current simulation time
19.     freq = neuron._frequency.x_   # Hz            The frequency of the oscillating activation signal
20.     amp = neuron._amplitude.x_    # unitless      Amplitude of the oscillating activation signal
21.     phasea = 0.
22.     phaseb = 180.
23. 
24.     # Calculate the oscillating activation signal
25.     adrive = amp*math.sin(2.*math.pi*freq*t + phasea*math.pi/180.)
26.     bdrive = amp*math.sin(2.*math.pi*freq*t + phaseb*math.pi/180.)
27.     neuron._MNLink1A_act.x_ = adrive
28.     neuron._MNLink1B_act.x_ = bdrive
29.     neuron._MNLink2A_act.x_ = bdrive
30.     neuron._MNLink2B_act.x_ = adrive
31.     neuron._MNLink3A_act.x_ = adrive
32.     neuron._MNLink3B_act.x_ = bdrive
33.     return

## Next time

We'll do one more example of a controller for this three-link pendulum. It will introduce the idea of center of mass control and introduce an external perturbation (i.e. we'll push the pendulum and see how it responds!)

The XML file is getting pretty large and cumbersome to work with. Fortunately we can use the embedded Python interpreter to programmatically build the model, run a simulation, and analyze the results, all within a single environment! This is nice because it streamlines the workflow of how most of us use biomechanical models. We're going to use this method to build our quadruped, but before we do that lets see how to do it with this three link model.