Skip to main content

Multibody Objects

The mBody object serves as a container for all the bodies, joints, forces, and sensors of a multibody system. It allows you to perform operations on the entire system at once - such as solving its equations of motion or animating results.

The mBody Constructor

To encapsulate bodies, joints, or other objects into a multibody system, type, for example,

>> mb = mBody(B, C, J)
mb = 
     _____             QuIRK: mBody object
    |\    \   /\\\\\\   2 bodies (28 states)
    | \____\o/  \\\\\\  1 joints (7 constraints)
    \ |    | \  /####/  0 sensors
     \|____|  \/####/   Equation of motion not yet solved.

Before you can solve, animate, or extract results from a set of bodies and joints, you must put them into an mBody container class. The mBody object includes functions that let you display and modify the entire system, as well as various analysis functions.

The mBody function takes as arguments a list of bodies, joints, forces, and sensors. You may list those objects in any order you like, followed by property/value pairs. The variable mb will then refer to the system as a whole. You can still get properties from and modify the individual objects inside the mBody object, but now you can use mBody functions as well.

mBody Options

Here is the full list of mBody constructor options. Give the constructor properties as strings, or use the shortcut abbreviations shown (they are not case-sensitive).

Damping Double N s/m Damping coefficient to apply to linear and angular velocities of all bodies in the system. Each body experiences a force -damping*v and a torque -damping*ω/radius.
Data Any user-defined data
Force Function handle Newtons Force function acting at all bodies' center of mass, as a function of position, attitude, velocity, angular velocity, and time: f(rqvωt)
Thresh Double m/s or N Threshold on zero-velocity, zero-force, and equilibrium detection in solutions
Torque Function handle N m Force function acting at all bodies' center of mass, as a function of position, attitude, velocity, angular velocity, and time: τ(rqvωt)
U Function handle J Potential energy of bodies in the system, as a function of either position and time U(r, t) or body object U(body)

You can obtain a list of all the mBody properties that you may later access of modify by entering the commands get(mBody) or set(mBody).

Back to top

Applying Body Forces and Torques

A multibody system will not move unless there are forces or torques acting upon it. QuIRK has four ways of specifying the forces and torques acting on bodies: individual-body forces, full-system forces, full-system potential energy functions, and inter-body force objects. When you solve an mBody, QuIRK calculates the total force on each body as the superposition of all four contributions.

Individual Body Forces and Torques

Each body has a force and torque property. These properties accept function handles that return the force vector and torque vector on that individual body's center of mass, as a function of the body's instantaneous state. That is, f(rqvωt) and τ(rqvωt). These functions may be simple Matlab function expressions, such as body1.force = @(r,q,v,w,t)(-0.5 * v), or may be a user-defined function that is much more complex.

Full System Forces and Torques

You can specify the force and torque on all body centers of mass in a multibody system in an analogous way, by setting an mBody's force and torque properties to f(rqvωt) or τ(rqvωt).

Potential Energy Functions

If you specify a potential energy property U on an mBody, then all bodies contained in that multibody system will experience a force f = -U, in addition to all explicitly defined forces. If the potential energy is independent of body properties except for position, then you should use the functional form U(r, t). However, if the potential energy does depend on body properties, the function U can accept as its argument a single body object. For example, setting

>> mb.U = @(B)( B.mass * 9.8 * B.pos(3) )

will produce the gravitational force -mg in the inertial z direction on all bodies in the multibody system, for g = 9.8 m/s.

If you plan to use the energy analysis functions after solving the mBody, this is the best way to express forces because it will make conservation of energy explicit.

Force Objects

Force objects that are members of the mBody define forces that act between bodies, depending on the state of each body. They produce equal and opposite forces and torques on the two bodies as given by the forceFn and torqueFn arguments to the force function.

Back to top

Solving mBodies

All the interesting stuff happens when you have defined a multibody system and want to solve its equations of motion. To find a solution for a multibody system, use the command

>> mb.solve(tspan)

where tspan is either a two-element vector containing the start and end times of the solution, or a vector of times. (This is the same format as the timespan input to Matlab's ode45.) The initial conditions for integration are the body states when you call solve, including initial velocity states - you can use mb.draw to check those states, and use commands like mb.settime to quickly jump to a state in a solved mBody (see "Modifying mBodies"). QuIRK will solve the multibody system's equations of motion and display a message when its solution is complete. As long as there are no hang-ups in the integration, it's as simple as that!

Once solved, the mBody object will contain its state history in the variable x and the time at each point in the solution in the variable t. If you display the multibody object in the command window, it will tell you the time range of the solution. After it is solved, an mBody object's state will be the final state in the solution history.

For example, consider the following sample commands from the file demo_Pendulum.m:

>> base = body([0 0 0.5], [0 0 0 1], 'shape', 'cube', ...
     'size', '2x2x1', 'color', 'm');
>> rod = body([0 2*sind(30) -2*cosd(30)], ...
     [sind(-30/2) 0 0 cosd(-30/2)], 'shape', 'cylinder', ...
     'mass', 0.1, 'size', '0.1x0.1x4', 'color', 'b');
>> bob = body([0 4*sind(30) -4*cosd(30)], ...
     [sind(-30/2) 0 0 cosd(-30/2)], 'shape', 'sphere', ...
     'size', '1x1x1', 'color', 'b');
>> grnd = joint(base, 'ground');
>> fixbob = joint(rod, bob, 'pt1', [0 0 -2], ...
     'pt2', [0 0 0], 'type', 'fix');
>> hinge = joint(base, rod, 'pt1', [0 0 -0.5], ...
     'pt2', [0 0 2], 'type', 'hinge', 'axis', [1 0 0]);
>> mb = mBody(base, rod, bob, grnd, fixbob, hinge, ...
     'U', @(b)( 9.8*b.mass*b.pos(3) ), 'damping', 1 ) 
mb = 
     _____             QuIRK: mBody object
    |\    \   /\\\\\\   3 bodies (42 states)
    | \____\o/  \\\\\\  3 joints (20 constraints)
    \ |    | \  /####/  0 sensors
     \|____|  \/####/   Equation of motion not yet solved.
>> mb.solve(0:0.1:10) 
-- QuIRK: Solving equations of motion...
              Solution complete in 2.0869 seconds.
ans = 
     _____             QuIRK: mBody object
    |\    \   /\\\\\\   3 bodies (42 states)
    | \____\o/  \\\\\\  3 joints (20 constraints)
    \ |    | \  /####/  0 sensors
     \|____|  \/####/   Equation of motion solved for t = 0:10.

Solver Options

You can provide the solve command with a set of properties denoting conditions at which the integration should halt, specify the Matlab integrator to use, and specify ode options. The four properties are:

Odeopts Structure Output of odeset ODE options structure
Sense Logical or Integer Stop simulation on sensor proximity
Solver String ODE function Choose ODE solver function
Stop String
Stop simulation when v ≈ 0, F ≈ 0, a weak test for collisions is true, or the specified number of realtime seconds elapses
Tidy Logical Perform extra state cleanup operations inside integrator

The 'sense' property will cause the mBody solution to halt if any sensor objects come within close proximity of one another (the key function of sensors). You can give an integer for 'sense'; doing so will cause solve to stop integration when at least the given number of sensors is in proximity. Sensors will only stop the integration when they come into proximity; sensors that start in proximity when integration starts will not cease integration. Note that solving with 'sense' enabled will decrease the value of RelTol in the integrator, which may slow integration. (If your sensors seem to pass through each other without triggering, you may need to decrease RelTol below 1e-8.)

Giving the property 'stop' the value 'zerovelocity' will cause the solution to stop when the velocities of all bodies in the system fall below the mBody's threshold value. Setting 'stop' to 'zeroforce' will cause the solution to stop if the forces on all bodies, in the direction of virtual displacements compatible with the constraints, fall below the threshold value. The 'collision' option will stop the solution when checkCollisions is true. Finally, a numeric value for 'stop' will cut off the integration after the number of seconds you specified (this option is good for situations when joints might lock up after some amount of time and you wish your code to continue). You can give the 'stop' property multiple times in the call to solve to enable both conditions.

Specify the Matlab integrator by putting 'solver' in the argument list, followed by the string name of the Matlab integrator function. The integrator must take the same inputs and provide the same outputs as, for instance, the default integrator ode45. For example,

>> mb.solve([0 5], 'stop', 'zeroforce', 'solver', 'ode23T');

finds a solution to the system's equations of motion from t = 0 to 5 s, using the Matlab integrator ode23T, and the solution will stop if the system experiences zero force in the direction compatible with the constraints.

Solve can also take as a property an ode options structure, as produced by the command odeset, which it will pass on to the integrator function. So, for example, to see a plot of the solution state as Matlab integrates, use

>> mb.solve(tspan, 'odeopts', odeset('OutputFcn', @odeplot));


When bodies move quickly relative to one another or get into certain configurations (e.g. near kinematic singularities), it is possible for numerical integration to generate slight violations of the joint constraints. This can cause the Matlab integrator to lock up. Set the 'tidy' property to true to instruct QuIRK to perform some extra clean-up operations each time the integrator calls the equation of motion. This may slow integration down, but produces more accurate and more robust results for longer integrations (greater than a few seconds with a few bodies).

Finding Equiibria

QuIRK can have Matlab stop its integration when the forces on a multibody system fall to zero (the system is at equilibrium). However, if we are looking to find that equilibrium state, we do not necessarily know how long to simulate for. The command

>> x0 = mb.findeq;

uses repeated calls to the mb.solve command to locate an equilibrium state even when we do not know in what time range the system is at equilibrium. It functions best when there is damping in the system. If you enter this command, QuIRK will print a progress bar as it searches for equilibrium and notify you when it finishes or reaches a pre-set integration time limit. When findeq finishes executing, it restores any previous state history in the mBody object.

>> x0 = mb.findeq;
     10 calls to mBody.solve.

The state x0 is in the same format as an mBody state vector. So, you can use mb.overwriteState(0, x0) to replace the mBody's state history with the equilibrium value, or simply use mb.draw('state', x0) to see what it looks like without overwriting any previous solutions.

Back to top

Modifying mBodies

You can modify many mBody properties by using the set/get or "." operator interface, such as the damping value of the multibody system. However, it is important for QuIRK to keep track of which bodies, joints, forces, etc are in which mBody object. So, if you want to add a new body, joint, or any other object to a multibody system, use the command

>> mb.add(newBody, newJoint, anotherNewBody);

You can supply the add command with a list of any QuIRK objects, except for other mBodies, in any order. Similarly, do not just clear components of an mBody from the workspace to erase them. Remove objects from an mBody with the command

>> mb.remove(oldJoint, oldBody);

Beware: these, and other, modifications to an mBody will clear any previous solution of its equations of motion.

After an mBody has been solved, it contains its full state history. If you want to set the mBody to its state at a particular time, use

>> mb.settime(timeIndex);

to assign each component object of the mBody to its state at the time given by the timeIndexth entry of the time vector stored in mb.t. If you do not specify any argument on settime, it defaults to the most recent time in the mBody's state history.

If you want to solve the multibody system again from the current state of all its subordinate bodies, you will have to clear the previous state history. Do this with the command

>> mb.setnow;

which will both clear the previous solution and set the mBody up for another call to solve. This is an important command because it tells the mBody to construct its new initial state vector by querying each component body for its position and velocity states. So, if you change the properties of any component bodies (see Modifying Bodies), or add or remove joints, and want to use those changes in subsequent simulations, you should call mb.setnow before starting the next solution. (Important note: Do not confuse setnow with settime, given no arguments. Settime just assigns each body to its state at the final time in the solution. Setnow will save the current time and clear the state history of an mBody.)

A quick shorthand is

>> mb.reset;

which acts like mb.settime(1) followed by mb.setnow.

Back to top

Drawing mBodies

You can display an mBody in much the same way as you would draw a body object. However, now there are a few extra options to play with. For example, you have already seen that on mBodies, the draw command supports an option called 'state' that lets you draw the system at a particular state, even if that state is not in its state history. The mBody draw function also lets you specify an index into a solved mBody's time vector, and will draw the mBody at that time index. In addition, you can specify the option 'filename' and draw will generate an output image depicting the multibody system. Ghost will work similarly on mBodies. See the draw help text for full options.

Besides draw and ghost, which work any time on an mBody and update "live" as you change body, joint, or mBody properties, you can view the trajectory of a solved multibody system.

First, the command

>> mb.animate

will produce a movie animation of the solved mBody in the active figure window. The movie will play in real time by default (unless the timesteps are too close together compared with the time it takes Matlab to draw everything and pause between frames). You can specify options on the animate command similar to the options on draw, such as lighting the figure or outputting an AVI movie. Refer to the animate help text for all the options. (Note: for best results when saving an AVI file, use an explicit, uniform time vector in mb.solve.)

Here is a sample movie of a double pendulum in gravity, with an initial velocity condition:

Second, the command

>> mb.trace

produces a translucent image of the mBody at intervals along the solution state history, with a more opaque image overlaid at the final state. This can be useful for seeing the extent of an mBody's motion, or seeing where the system trajectory spends more or less time.

traced mbody
Back to top

Analysis Functions

The mBody class contains built-in functions to access and plot the momentum and energy of a system after it has been solved. The commands

>> [K U] =;


>> [p H] = mb.momentum;

return the instantaneous kinetic and potential energy or linear and angular momentum vectors of the system, respectively, in its current state. (These two commands also work before an mBody is solved.)

You can create a plot of the total system energy as a function of time, with the command

>> [K U] = mb.plotEnergy;
traced mbody

This command also returns the time histories of the kinetic and potential energy of the system in the variables K and U. They have the same number of elements as the vector mb.t. (This command is where specifying forces on an mBody as a potential energy gradient pays off). Notice in the example plot, from the pendulum system, that the mBody's damping value always makes the total energy decrease. The change in energy might be something interesting to look at, so QuIRK also includes a command to plot and save the total external power applied to the system:

>> P = mb.plotPower;

There is a similar command,

>> [p H] = mb.plotMomentum;

which plots a time history of the total magnitude of the linear and angular momentum of each body in a solved mBody.

Back to top



  • Joseph Shoer

Download QuIRK

Requires MATLAB 2009 or later