Extending the pendulum to three links

From neuromechanicwiki
Jump to: navigation, search

In the Building your first Neuromechanic Model tutorial we constructed a pendulum and simulated its periodic motion. The model had two bodies, a fixed base and the pendulum. We will now extend the model to a three-link pendulum and use python or matlab to do some simple analysis of its motion.

Extending the model to three links

Recall that our model looks like this:

<NeuromechanicFile>
   <Bodies>
      <RigidBody name="base">
         <FrameLocation Parent="ground">0 0 0</FrameLocation>
         <DegreeOfFreedom Name="base1tx" Type="translation">
            <Axis>1 0 0</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>true</Locked>
         </DegreeOfFreedom>
         <DegreeOfFreedom Name="base1ty" Type="translation">
            <Axis>0 1 0</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>true</Locked>
         </DegreeOfFreedom>
         <DegreeOfFreedom Name="base1tz" Type="translation">
            <Axis>0 0 1</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>true</Locked>
         </DegreeOfFreedom>
         <DegreeOfFreedom Name="base1rx" Type="rotation">
            <Axis>1 0 0</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>true</Locked>
         </DegreeOfFreedom>
         <DegreeOfFreedom Name="base1ry" Type="rotation">
            <Axis>0 1 0</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>true</Locked>
         </DegreeOfFreedom>
         <DegreeOfFreedom Name="base1rz" Type="rotation">
            <Axis>0 0 1</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>true</Locked>
         </DegreeOfFreedom>
         <Mass>0.00418</Mass>
         <Inertia>1.6755e-007 1.6755e-007 1.6755e-007 0 0 0</Inertia>
         <CenterOfMass>0 0 0</CenterOfMass>
      </RigidBody>
      <RigidBody name="pendulum">
         <FrameLocation Parent="base">0 0 0</FrameLocation>
         <DegreeOfFreedom Name="pendulum1z" Type="rotation">
            <Axis>0 0 1</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>false</Locked>
         </DegreeOfFreedom>
         <Mass>0.1676</Mass>
         <Inertia>3.4851e-004 2.6808e-005 3.4851e-004 0 0 0</Inertia>
         <CenterOfMass>0 0.12 0</CenterOfMass>
      </RigidBody>
   </Bodies>
</NeuromechanicFile>

Since we are going to add a few more links to the pendulum I'm going to change the name of the second rigid body from "pendulum" to "link"

      <RigidBody name="link1">
         ...
         <DegreeOfFreedom Name="link1z" Type="rotation">
         ...
      </RigidBody>

Now all we do, is duplicate this body with a few changes

  • The new rigid body needs a unique name, lets call it "link2"
  • The new rigid body degree of freedom needs a unique name, lets call it "link2z"
  • The new rigid body will be attached to the top of "link1". We'll put it at a distance of 24 cm from link1's origin in the y-direction.
      <RigidBody name="link2">
         <FrameLocation Parent="link1">0 0.24 0</FrameLocation>
         <DegreeOfFreedom Name="link2z" Type="rotation">
            <Axis>0 0 1</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>false</Locked>
         </DegreeOfFreedom>
         <Mass>0.1676</Mass>
         <Inertia>3.4851e-004 2.6808e-005 3.4851e-004 0 0 0</Inertia>
         <CenterOfMass>0 0.12 0</CenterOfMass>
      </RigidBody>

Lets repeat this one more time with a link named "link3" attached to "link2". You should get this:

<NeuromechanicFile>
   <Bodies>
      <RigidBody name="base">
         <FrameLocation Parent="ground">0 0 0</FrameLocation>
         <DegreeOfFreedom Name="base1tx" Type="translation">
            <Axis>1 0 0</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>true</Locked>
         </DegreeOfFreedom>
         <DegreeOfFreedom Name="base1ty" Type="translation">
            <Axis>0 1 0</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>true</Locked>
         </DegreeOfFreedom>
         <DegreeOfFreedom Name="base1tz" Type="translation">
            <Axis>0 0 1</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>true</Locked>
         </DegreeOfFreedom>
         <DegreeOfFreedom Name="base1rx" Type="rotation">
            <Axis>1 0 0</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>true</Locked>
         </DegreeOfFreedom>
         <DegreeOfFreedom Name="base1ry" Type="rotation">
            <Axis>0 1 0</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>true</Locked>
         </DegreeOfFreedom>
         <DegreeOfFreedom Name="base1rz" Type="rotation">
            <Axis>0 0 1</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>true</Locked>
         </DegreeOfFreedom>
         <Mass>0.00418</Mass>
         <Inertia>1.6755e-007 1.6755e-007 1.6755e-007 0 0 0</Inertia>
         <CenterOfMass>0 0 0</CenterOfMass>
      </RigidBody>
      <RigidBody name="link1">
         <FrameLocation Parent="base">0 0 0</FrameLocation>
         <DegreeOfFreedom Name="link1z" Type="rotation">
            <Axis>0 0 1</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>false</Locked>
         </DegreeOfFreedom>
         <Mass>0.1676</Mass>
         <Inertia>3.4851e-004 2.6808e-005 3.4851e-004 0 0 0</Inertia>
         <CenterOfMass>0 0.12 0</CenterOfMass>
      </RigidBody>
      <RigidBody name="link2">
         <FrameLocation Parent="link1">0 0.24 0</FrameLocation>
         <DegreeOfFreedom Name="link2z" Type="rotation">
            <Axis>0 0 1</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>false</Locked>
         </DegreeOfFreedom>
         <Mass>0.1676</Mass>
         <Inertia>3.4851e-004 2.6808e-005 3.4851e-004 0 0 0</Inertia>
         <CenterOfMass>0 0.12 0</CenterOfMass>
      </RigidBody>
      <RigidBody name="link3">
         <FrameLocation Parent="link2">0 0.24 0</FrameLocation>
         <DegreeOfFreedom Name="link3z" Type="rotation">
            <Axis>0 0 1</Axis>
            <State Type="Kinematic">0 0</State>
            <Locked>false</Locked>
         </DegreeOfFreedom>
         <Mass>0.1676</Mass>
         <Inertia>3.4851e-004 2.6808e-005 3.4851e-004 0 0 0</Inertia>
         <CenterOfMass>0 0.12 0</CenterOfMass>
      </RigidBody>
   </Bodies>   
</NeuromechanicFile>

I named my file "threelink.nmcb". Save yours and load it into Neuromechanic

Simulating motion of a 3-link pendulum

Before we begin the simulation lets change a few things. First we need to add gravity, switch to Source view and find the <Gravity> element (it will have been automatically created and initizialized to zero. Set it as below.

<NeuromechanicFile>
   <Bodies>
      ...
   </Bodies>   
   <Environment>
      <Gravity>0 -9.80665 0</Gravity>
      ...
   </Environment>   
</NeuromechanicFile>

In the Model Tree click on bodies to open the degree of freedom editor.

  • Change the posture of the three links so that the model is a configuration that will move under the influence of gravity (not straight up or straight down).
    • I used 7, 17, and 23 degrees respectively for link1z, link2z, and link3z

In the Model Tree click on parameters to open the simulation editor.

  • Change the End Time of the simulation to 3.
  • Change the Reporting Interval to 0.001
  • Select the Save Output checkbox
  • Select output variables ModelCOM and ModelKineticEnergy.

Press Simulate in the simulation toolbar.

Analysis of motion

It's nice to look at simulation but one of the most powerful things about models is that we can get information out of them that is very difficult to obtain from experiments. Lets demonstrate that by looking at the energy of the pendulum.

xml2matlab

We'll do this in matlab (or Octave). You will need a file (xml2matlab) available here which takes as an argument the file name to be parsed (make sure it has the appropriate path or is in the current working directory) and returns a cell structure of the contents of the xml file. xml2matlab is a general xml parser that works for both Matlab and Octave. However, it does treat the elements a bit differently than many parsers. If the content of the element is numerical (scalar or array) then it will automatically convert it from a character string to an array. It also saves the order in which each element was read.

The first time you use xml2matlab you might want to examine the structure to make sure you understand how the data is arranged. Here is function that can help you do that if you like. Step through the code line by line to explore the entire structure.

% Convert xml to a matlab cell structure
xml{1} = xml2matlab('myfile.nmcb');
printxmlstruct(xml,'xml','')

where

function printxmlstruct(A,name,buffer)
    if iscell(A); % This is an element with children
        n = length(A);
        for ii = 1:n
            if isstruct(A{ii})
                disp([buffer '<',name,'>' num2str(ii) ' of ' num2str(n)])
                FN = fieldnames(A{ii});
                for jj = 1:length(FN)
                    printxmlstruct(A{ii}.(FN{jj}),FN{jj},['  ' buffer])
                end
            elseif strcmpi(name,'Comments__')
 
                disp([buffer 'Comment ' num2str(ii) ' of ' num2str(n) ':' A{ii}])
            else                warning('Something failed')
            end
        end
    elseif isstruct(A) %This is an attribute
        disp([buffer,name,'='])
        FN = fieldnames(A);
        for jj = 1:length(FN)
            printxmlstruct(A.(FN{jj}),FN{jj},['  ' buffer])
        end
    else %This is the text node
        if ischar(A)
             disp([buffer name '=' A])
        else
            disp([buffer name '=' mat2str(A)])
        end
    end

You can use any xml parser you like but you'll need to change the syntax of the code below to match your parser.

Mechanical energy analysis

Threelinkenergy.png

So what happens with the energy in the pendulum? If you remember your classical mechanics, you'll remember that the total mechanical energy of a system should sum to a constant at all points in time (conservation of energy) IF there are no non-conservative forces acting on the system. The conservative forces can be accounted for with potential energy. In our case the only conservative forces are due to gravity and we can use <ModelCOM> to track the potential gravitational energy during the simulation. We also output the mechanical kinetic energy in <ModelKineticEnergy>.

nm = xml2matlab('threelink.nmco');                  % Convert xml to a matlab cell structure
bod = nm.NeuromechanicFile{1}.Bodies{1}.RigidBody;  % Get <RigidBody> elements
mass = 0;                                           % Initialize mass to zero
for ii = 1:length(bod)
    mass = mass+bod{ii}.Mass{1}.Value;              % Add total mass from each body
end;
dyn = nm.NeuromechanicFile{1}.Dynamics;             % Get <Dynamic> elements
g = nm.NeuromechanicFile{1}.Environment{1}.Gravity{1}.Value';   % Get gravity vector
n = length(dyn);                                    % Get the number of data points
t = zeros(n,1);                                     % Initialize time
KE = zeros(n,1);                                    % Initialize kinetic energy
COM = zeros(n,3);                                   % Initialize center of mass position
PE = zeros(n,1);                                    % Initialize potential energy
for ii = 1:n
    t(ii,1) = dyn{ii}.Time{1}.Value;                % Store time
    KE(ii,1) = dyn{ii}.ModelKineticEnergy{1}.Value; % Store kinetic energy
    COM(ii,:) = dyn{ii}.ModelCOM{1}.Value(1:3);     % Store center of mass position
    PE(ii,1) = -(COM(ii,:)-COM(1,:))*g*mass;        % Calculate potential energy w.r.t initial CoM Position
end;
plot(t,[KE PE KE+PE]);                              % Look at energy exchange
xlabel('Time (sec)')
ylabel('Energy (J)')
legend({'Kinetic Energy' 'Potential Energy' 'Total Energy'})

If you did this like I did you should see the three traces above. Energy is conserved!

In the next tutorial we'll add muscles to