What’s better suited for object-oriented programming (OOP) than building a simulation of a physical system? Each physical object (e.g. mass, spring etc.) can be implemented as an object and then we should be able to easily build and modify large models.
However, there are a lot of challenges to get there, including:
- How to work with classes and objects at all (if you’re like me new to OOP)?
- How to build a working object structure that adopts to the requirements of a typical ODE solver (which usually requires access to all state variables in the so-called state-space representation)
- How to design classes for different object types in a way that allows them interact with each other to calculate and transmit forces and accelerations.
Getting there turned out to be quite interesting and I finally got a simple object-oriented multibody model to simulate in Matlab and GNU Octave. This post is a wrapup of all the stuff I learned to get there.
Ressources: Getting started with object-oriented programming in Matlab or GNU Octave
Normally, I try to keep my blog posts as much as focused towards GNU Octave as possible to make them acessible for everyone (while I have access to Matlab both at my university and work desks, I do not have a licence for my home computer and I do not want to exclude all those readers who are in the same situation). I very much enjoy the fact that Octave is available as free software. However, here I need to praise Matlab for the amazing documentation they provide along with their products, which turned out to be the by far best ressource to learn about all the object-oriented stuff: See the Matlab online help system or get the PDF with the complete guide to object-oriented programming (800 sites of detailed and easy-to-navigate information, might require a (free) Mathworks account that’s also needed for the online tutorials and Matlab Central). I mainly read parts of the PDF to get started and found that it easily answered all the questions I had. Most of it also works in GNU Octave but you need to figure that out on your own.
Additional great ressources are the introductions to object-oriented programming in Matlab given here and here. Depending on your focus, you might prefer them over the official Matlab documentation (especially if you don’t want to dig through 800 pages of pdf material).
Additionally, when you are working with objects in Matlab or GNU Octave, you might want to know these commands:
- use whos to see information on all variables (including your objects) in the current workspace
- use clear all to clear all objects and the currently loaded class definitions (use this if you changed some of the code in one of your classes to make sure you’re not working with some objects built from the old class).
- use isa(object, ‘classname’) to check if <object> is part of a class. Example: isa(mass_1, ‘handle’) returns true if mass_1 is a child of the handle superclass.
If you are completely new to object-oriented programming (like me), you should first learn about these concepts:
- class vs object (a class defines the behavior for all objects of this class, e.g. mass_1d is a class that defines the general behavior for individual objects m1, m2, … we’ve built on it), see “classes and objects” here or page 1-3 in the pdf.
- handle objects (the handle points to the object, which is stored somewhere else in your computer’s memory – when you copy a handle variable, both handles point to the same object), see pdf chapter 7 or here.
- class inheritance (classes can be built upon other classes and inherit all their features), see pdf chapter 12 or here.
Simulation model structure
My approach is this: To create a multibody simulation with an underlying object structure in Matlab, we need three things:
- Multibody objects, defined as classes in the corresponding Matlab OOP syntax.
- An object stack that contains all the objects our multibody system is made of.
- A wrapper function that communicates between the numerical integrator (e.g. ode45) and the object stack.
Our object stack might be a global struct variable containing all objects and maybe additional information. The wrapper function needs to supply access to the object stack for the solver as a state-space representation function – something like this:
yp = function(t,y) % yp = dy/dt, y = state vector, t = time % 1) set state y in object stack % 2) calculate derivative yp for each object, somehow like: % 2a) update all objects % 2b) calculate yp for each object % 3) build yp vector based on objects
This sounds very simple but it took me a lot of thinking to get there. Good questions are: How will my objects know what to do? Let’s say I have a spring connected to 2 masses. At 2a in the list above, the spring needs to a) get the current state from the mass it’s connected to, b) calculate the spring force based on the positions and the stored spring rate and c) transmit the calculated force back to to objects it’s connected to. How will the spring know which object’s it’s connected to? I guess we can store the connection to those as handles, or maybe indexes, but we’ll need to discuss this later on…
My current implementation
Most of the choices I have made are already hinted at above. My implementation currently features:
- A global struct variable s that holds all the objects for my multibody model (holds both the object handles and the information on which part of the state vector corresponds to each object).
- As a standard interface for all objects: An (abstract and handle) parent class MBS_object that contains state and parameter variables and defines all the standard functions the programming structure expects: set_state, calc_int_state and calc_yp
- Example classes for my multibody system: mass_1d and spring_1d classes
- A wrapper function MBS_ODE that accepts the state-vector y from the ODE solver, updates the objects and outputs the current yp vector
The program structure is given below:
(note that I am basically clueless on how to actually draw UML diagrams or whatever. This seems to represent what my code contains)
The spring class also stores handles to the two masses it is connected to. During the update call, it grabs the position from both masses, calculates the spring force and adds the force back to the mass objects:
Note that currently, all object parameters are set to global – this allows the spring object to directly access the mass object’s state in line 30 above. With larger models, setting the state vector to private and adding a get_state function might be a good idea.
The ODE_MBS function basically calls the three functions defined in the MBS_object superclass:
To get a convenient connection between state vector and internal state, the vector i that also sits in the s structure is used (it contains the indexes of the state variables in the solver’s state vector). This allows easy passing of values between ODE solver and objects:
- set_state(object, y( indexes ) ) pushes the state vector into an object.
- yp( indexes ) = calc_yp( object ) pulls the yp part for an object back into the yp vector for the ODE solver.
Results and Discussion
The create_simulate_MBS-script sets up a few objects and simulates them. Currently, there is no automatic creation of the i-parts in the struct (I’ve just manually added them). However, creating objects and simulating them is quite simple and works with a few lines of code:
While building the object-oriented structure itself was a lot (!) more effort than with creating a simple state-space representation for the same system, building different models out of the defined classes is very simple.
The results also look as expected. Here are the simulation results for a simple 3-mass system:
Next steps might be:
- Adding more objects
- Improving the creation of new objects, for example with an automatic creation of the i vector in the objects stack
- Automatic creation of a model based on a text-file input
- Porting the whole software to C++ to get better portability and highly reduced simulation times.