Python

Python is a general-purpose and interpreted programming language, first developed by Guido van Rossum and released in 1991. In order to make it freely accessible to everyone, it is developed under an open-source license. Python Software Foundation, which is responsible for administering Python’s license, was formed “to promote, protect, and advance the Python programming language, and to support and facilitate the growth of a diverse and international community of Python programmers” 28. Due to its high-level and open-source nature, large community of developers in Python, an abundance of packages developed for Python 29, and easy integration with other programming languages, Python has gained popularity in various applications and disciplines and has turned into an industry standard. Consequently, many valuable resources have been developed for getting started with Python. The website 28 provides a rich repository of documentation about Python installation, learning materials, news and events. In addition to many books written about Python 30, 31, several MOOCs (Massive Open Online Courses) also exist to provide familiarity with the fundamentals, advanced topics, and applications of Python 32, 33.

In the past decade, the Python community has developed numerous packages for mathematics, science, and engineering applications 34, 35. The following is a list of some of these libraries

• SciPy: A comprehensive library with sub-modules for various mathematical and scientific operations such as Fast Fourier Transform, interpolation, numerical integration, linear algebra, file input/output, optimization and curve-fitting, statistics, and signal processing.

• NumPy: Offers computational capabilities similar to Matlab which also overlaps with some of SciPy’s functionalities. It is best suited for fast array operations. Matlab users who wish to learn more about NumPy can refer to 36 for a detailed comparison between Matlab and NumPy.

• Matplotlib: Enables high-quality 2D plotting. It can be used to generate plots, histograms, power spectra, bar charts, error charts, scatterplots, etc.

• SymPy: Used for symbolic programming. It is incorporated into a wide range of other libraries such as symbolic statistical modeling, multibody dynamics, linear circuit analysis, etc.

• IPython: Provides an interactive shell similar to Matlab’s workspace which is ideal for troubleshooting, interactive data visualization, use of GUI toolkits, and parallel computing.

• pandas and Scikit-learn: Offer an extensive suite of tools for predictive data analytics, modeling, and machine learning.

• Python Control Systems Library: Provides functions and classes to systematically analyze and design feedback control systems along with a Matlab compatibility module for users more familiar with Matlab 37.

As mentioned earlier, Python has received an increased attention in the past few years, mainly due to the fact that while it is a general-purpose, high-level, and interactive programming language, it is also equipped with powerful features and add-on libraries which place it at the same high level as specialized software such as Matlab. To catch up to its fast pace of growth and popularity and to familiarize the future generation, Python is increasingly being introduced in higher education curricula. In the past decade, majority of top CS departments in United States have switched to Python in their introductory computer programming courses 38. According to a comparative study in 39, among Python, C, and Matlab as the programming language of choice for engineering education, Matlab is determined to be more suitable compared to C, however, Python is selected as the best option due to its clarity and functionality. A similar conclusion about Python being the best choice for engineering education compared to Matlab and C is achieved in 2. Authors in40 discuss the motivation and the process of a transition from C to Python in the first computer programming course for Mechanical Engineering students and report a reduced failure rate among the students, among other observations. In the context of MRE education, Python has mainly been used as the programming language for open-source hardware such as Raspberry Pi 41, 42, which is typically used for teaching control engineering. Python has also been used in43 to develop internet tools to teach the fundamentals of control.

Java

While Java is not the first programming language environment that comes to mind in an MRE context, it has many properties that make it a very useful tool. The first of these is the language itself. Java is fully object-oriented, has a relatively clean syntax, executes quickly, and operates with no changes (usually) across different computers and operating systems. It can be applied to mathematically based problems (for example, dynamic simulation) as well as used for direct control of physical systems.

Java compiles to an abstract machine code that executes via a Java Virtual Machine (JVM). Thus, a compiled Java program can execute on any system for which a JVM exists. This guarantees that basic properties such as how integer and real numbers are stored, precision associated with different data types, etc., will not change from one environment to another.

Java is open source. The Java Development Kit (JDK) including, compiler, linker, etc. is available from 44. Commercially licensed versions are also available from Oracle. The two most popular Integrated Development Environments (IDE) are also open source, Netbeans 45 and Eclipse 46.

While the Java language itself provides very effective ways to organize complex problems, solving them often requires additional numerical math tools. For a wide range of problems, the appropriate tools are available through the open source Apache Commons Math Library 47 or its successor, the Hipparchus Math Library 48. There are also other smaller, more specialized libraries that can be found online.

Java is used as the computing language for the Advanced Placement Computer Science course and is also widely taught in lower division undergraduate courses in computer science as well as other programs. Because of this, there is a large amount of material on learning Java – so much that any listing here would barely scratch the surface. The material consists of books, online tutorials, videos, forums, etc., much of it freely available. In general, this material focuses on learning Java without going very deeply into the subject matter used for examples. For material that is more closely focused on the subject matter of this contribution, a well-documented course on using Java for simulation of physical systems can be found at 49. A more extensive discourse on the use of Java in simulation has been produced under the auspices of the Open Source Physics project 50.

Modelica

“The Modelica Language is a non-proprietary, object-oriented, equation-based language to conveniently model complex physical systems containing, e.g., mechanical, electrical, electronic, hydraulic, thermal, control, electric power or process-oriented subcomponents” 51.

Modelica is of particular interest to the mechatronics/robotics community because it provides dynamic simulation of three-dimensional and constrained mechanical systems and energetically correct interactions between the mechanical systems and electrical systems, hydraulic, pneumatic, etc. It does this through an underlying mathematical base that is capable of solving differential-algebraic equations (DAEs) and an underlying computing syntax based on equations rather than algorithmic “store-into” statements. Most users do not interact with Modelica at that level, however, but rather make use of an extensive set of libraries for various physical domains combined with a graphical user interface.

Systems with constraints are known as “acausal”. For a system defined by a set of N first-order, nonlinear differential equations, constraints add algebraic equations to the system equations. Acausal systems have fewer than N independent initial conditions because the algebraic constraints reduce the effective order of the system. In mechatronics, almost all mechanical systems of any interest have constraints which is one reason Modelica is so valuable.

The Modelica language is defined by an open source specification maintained by the Modelica Association [modelica.org]. The actual implementation is left to any organization, including both commercial and open source entities. The most important open-source implementation is OpenModelica 52.

Modelica has been heavily used in certain industries, but less so in college-level teaching and, probably, not at all at the high school level. Thus, there is not nearly as much instructional material available as there is for much more widely used languages such as Python or Java. There is, however, some useful material to cite. The Open Source Modelica Consortium53 has a short, introductory course, as well as a longer, more in-depth course at 54. “Modelica by Example” is an in-depth, online book on Modelica 55. It is also available for sale (on a “pay what you can” basis) in various eBook formats. A set of introductory slides that uses many examples from “Modelica by Example” is available at 56. Because Modelica is explicitly aimed at simulation of physical systems, any instructional material on Modelica is, by definition, also about simulation of physical systems.

Gazebo

Since its debut in 2007, ROS (Robot Operating System) has grown into the most prominent open-source middleware framework for robotics research and development in industry and academia. Likewise, the Gazebo Simulator has become increasingly popular, thanks in part to its robust integration with the ROS framework 57. Gazebo is an open-source 3D dynamic simulator that provides a framework for realistic simulation of multiple robots dynamically interacting with each other in complex environments and scenarios 58, 59. Gazebo comprises a suite of features which include high-performance physics engines (e.g. Open Dynamics Engine, Bullet, etc.), high-fidelity graphical rendering of environments and ability to generate realistic sensory data from the environment 58. These features allow users to create very realistic simulated environments and scenarios to learn to operate and test their robot designs and control algorithms without the risk of damaging their physical robots. Users can either leverage the rich library of existing robot models available with the software (ranging from mobile aerial and ground robots to state-of-the-art humanoid robots) or design their own robot models from scratch using the Unified robot description format (URDF). As mentioned earlier, Gazebo will be used in the second part of this work to simulate the closed-loop performance of a robot manipulator.

In the context of MRE education, Gazebo and ROS have been used in several ways to help students learn core robotics concepts. A two-part introductory robotics course was presented in60 using both Gazebo and ROS. A graduate course on mobile robotics and robotic manipulation in the context of Industry 4.0 was developed in61 using both software. In 62, Gazebo and ROS were both applied in an advanced course on humanoid robotics. In robotics education, open-source simulation tools are particularly relevant because they address the challenge of hardware affordability, making it possible to design labs partially or entirely in simulation 61, 63. The integration of these open-source software into the curriculum requires that the students are provided with opportunities to learn how to use them. For instance, in 61, the lab included preparatory sessions where the students complete comprehensive tutorials on the software available online 64.

GNU Octave

GNU Octave 65 is a high-level programing language for numerical computations with built-in plotting and visualization tools. Octave is an open-source software with syntax largely comparable to and compatible with MATLAB. It can be freely redistributed and/or modified under the GNU General Public License. It runs on multiple platforms including GNU/Linux, macOS, BSD, and Windows. Octave was primarily developed by John W. Eaton and was conceived in 1988 as companion software for a chemical reactor design textbook. It is an interpreted programming language designed to solve a variety of engineering problems involving the numerical solution of linear and nonlinear problems and ordinary differential equation problems. It also contains an extensive collection of packages that enhances its core functionality by introducing field specific features. Some of the packages relevant to MRE systems include the Control, Image Acquisition, Image, Optim, and the Arduino package which are used for control design and analysis, image acquisition and processing, optimization, and Arduino control, respectively. Octave also includes both a graphical user and a traditional command-line interface that provides greater flexibility for novice and advance users. It can also be used to write non-interactive user-defined programs and can call C, C++, and Fortran code within Octave using Oct-Files or using MATLAB-compatible Mex-Files. For new users and developers Octave also provides an active Wiki-page containing installation, development, and documentation resources, in addition to the built-in Octave documentation 66.

In engineering education, GNU Octave has been used as an alternative platform to provide a computational environment to process and visualize data. Its appeal is linked to the fact that it is a free and open-source licensed product that is highly compatible with the MATLAB programming language. Universities have utilized Octave to create line following and robot simulators and to process data from mobile robot projects 67, 68, 69. It has also been combined with open source hardware such as the Raspberry Pi to perform image processing tasks in robot arm sorting applications 70. One of the main challenges with Octave that educators have pointed out is the lack of repository libraries like MATLAB’s mature toolboxes 68. Nevertheless, Octave has emerged as a promising alternative tool to MATLAB without compromising performance and productivity 70.

Motor Case Study

A DC motor is one of the most commonly used actuators for generating rotational and translational motions in MRE systems. Furthermore, DC motors are perfect examples of an MRE system due to their multidisciplinary nature. The DC motor dynamics are governed by linear differential equations, and therefore, familiarity with numerical simulation of DC motor dynamics can pave the way to study other MRE systems that the students and professionals may encounter in practice.

The circuit diagram below is a typical schematic of an armature-controlled DC motor.

In this DC motor, the armature voltage v, is the input to the system whereas the rotational displacement of the motor shaft θ, the angular velocity of the motor shaft ω = dθ/dt, and the armature current i are assumed to be the system outputs. Using Kirchhoff’s and Newton’s laws, the dynamics of this system can be described by the following set of differential equations:

where J is the moment of inertia of the motor shaft, b is the viscous friction coefficient of the motor shaft, K_t is the armature constant, L is the armature electrical inductance, R is the armature electrical resistance, and K_b is the back-emf constant. The values of these parameters used here are obtained from the datasheet information and model identification tests conducted on a brushed servomotor, C40-E-500FE, from Advanced Motion Controls 72. These parameters are reported in Table 3.

To conform to common numerical simulation algorithms required to solve the DC motor dynamics in Equation 1 , the equations first need to be converted to a set of first-order ODEs as shown below:

These equations can also be written in a state-space form:

where x(t) = y(t) = [θ(t), ω(t), i(t)]^T, u(t) = v(t), and the matrices A, B, C, and D are:

Python

In order to obtain the step response of the DC motor, the governing ODEs can be directly solved in Python using solve_ivp command from the scipy.integrate library 74. The other library required for this simulation includes numpy (imported as np). The solve_ivp command syntax for solving the ODEs in Equation 2 or Equation 3 is:

sol = solve_ivp(model, (t0, tf), x0, dense_output=True, vectorized=True, args = (v, ))  
x = sol.sol(t)

where model is a user-defined function. For the ODEs inEquation 2 , this function can be written as:

def model(t, x, u):    
    theta, omega, i = x  
    dxdt = [omega, -(b/J)*omega+(Kt/J)*i, -(Kb/L)*omega-(R/L)*i+(1/L)*u]       
    return dxdt 

whereas for the ODEs in Equation 3 , the function is:

def model(t, x, u):
    A = [[0, 1, 0],[0, -(b/J), (Kt/J)],[0, -(Kb/L), -(R/L)]]  
    B = [[0], [0], [1/L]]  
    dxdt = np.dot(A,x)+ np.dot(B, u)  
    return dxdt

In Python, the model parameters can be specified after the definition of the model function. The simulation time is defined using the arange command as 0 ≤ t ≤ 0.5 with a step size of 0.001 s. Finally, the input voltage v is chosen as the nominal motor voltage, i.e. 24 V, and the initial conditions vector, x0, is defined as a 1-dimensional list, x0 = [0, 0, 0].

Python Control Systems Library 37 can provide alternative methods for solving for the step response of the DC motor dynamics. The best way to install this library is through using a pip installer, more information can be found in the official documentation 37. After the installation, this library needs to be imported to the Python workspace (imported as ctrl here.) The step response can then be obtained using the step_response command from the Python Control Systems library:

t, x = ctrl.step_response(v*ssModel, t)

where ssModel, the state-space representation of the system described in Equation 3, is defined as:

 ssModel = ctrl.ss(A,B,C,D)  

While the classes and commands in Python Control Systems Library can be used for simulating the step response, Matlab compatibility module, also included in this library, can provide another convenient method for model simulation, especially for those more familiar with Matlab. Using this module, the step response can be obtained as:

 ssModel = matlab.ss(A,B,C,D)  
x, t = matlab.step(v*ssModel, t) 

where matlab is an alias used for importing the control.matlab module. Note that the state vector, x, obtained using Matlab Compatibility Module has a dimension of n × 3 whereas the previous methods generate a state vector with a dimension of 3 × n, where n is the length of the time vector. Finally, the library matplotlib.pyplot can be used for plotting system states and plot customizations. Figure 3 shows the step response of the DC motor, plotted using matplotlib.pyplot library in Python. Although slightly different compared to Matlab, this library provides a wide range of programmatic customization tools for plotting and modifying the appearance of the plots and is easy to implement 75. The generated figures, however, do not have a lot of the interactive functionalities such as setting cursors, graphical modification of the plots, etc, that Matlab provides.

Other system descriptions such as transfer functions and frequency response data models can also be defined in Python Control Systems Library using tf and frd commands, respectively. Other features of Python Control Systems Library which can be very helpful in the context of Mechatronics and Robotics Engineering can be found within the online documentation 37. Examples of some of these capabilities are given in the next section alongside similar features in Octave Control package.

GNU Octave

To obtain the step response of the DC motor to a constant armature voltage of 24 volts, two approaches are used in this section. The first approach uses Octave’s built-in function ode45 while the second approach implements the step function in Octave’s Control package to obtain the time response of the governing equations of the motor. For this example, it is assumed that the motor starts from rest.

For the first approach, the ode45 command syntax for solving the ODEs in Equation 2 orEquation 3 is:

 [~,X] = ode45(@(t,x)motor.model(t,x,v),time,x0,options);

where time is the time vector for the simulation, x0 is the initial condition for the system, options are the solver options, and model is an instance method from the user-defined DCMotor object named motor which implements the ODEs in Equation 2 as shown below:

 function xdot = model(obj,t,x,u)   
    theta = x(1); omega = x(2); i = x(3);   
    v = u; % armature voltage           
    xdot = [omega, (obj.Kt*i-obj.b*omega)/obj.J, ...   
            (v-obj.Kb*omega-obj.R*i)/obj.L];  
end   

The obj refers to the instance of the DCMotor class. To implement the state-space representation as in Equation 3 , the model method is defined as:

function [xdot] = model(obj,t,x,u)    
    A = [0 1 0; 0 -obj.b/obj.J   obj.kt/obj.J;0 -obj.kb/obj.L   -obj.R/obj.L];   
    B = [0;0;1/obj.L];  
    C = [1 0 0;0 1 0;0 0 1]; D = 0;         
    xdot = A*x+B*u;  
end  

For the second approach, Octave’s Control package 76 needs to be installed and then loaded. Once the package is available, the motor step response can be simulated using the following syntax:

ssModel = ss(A,B,C,D);  
[~,~,X] = step(ssModel*u,time);  

where the ss function creates an LTI state-space representation of the motor and the step function determines the step response of the LTI system. The simulated state vector is then stored in the variable X. Plotting these state variables results in a figure identical to Figure 3 .

Similar to Python’s Control Systems Library, Octave’s Control Package provides other system descriptions for defining transfer functions and frequency response data models. It also features other useful control system design and analysis tools. A brief summary of some of the available control-related commands in both Python and Octave are shown in Table 4.

Modelica — Equation Mode

Using Modelica in its equation mode comes very close to just transcribing the dynamic equations and clicking “GO”. Below is the content of the “equation” section of the Modelica model:

der(theta) = omega;  
der(omega) = ((-b * omega) + Kt * i) / J;  
der(i) = ((-Kb * omega) - R * i + v) / L;  

where der is derivative-with-respect-to-time. Note that these are true equations, not computing statements. They can be written in any order and, within each equation, the terms can be rearranged as long as the original algebraic meaning is maintained.

These listings are from OpenModelica’s Connection Editor. OpenModelica was used for this case study and all Modelica examples to follow 52.

As with any problem definition, the parameters and initial conditions must also be specified. This is also straightforward in Modelica. The parameters are given by:

parameter SI.Resistance R = 2.62 " Ohm";  
parameter SI.Voltage v = 24.0 " v";  
parameter SI.ElectricalTorqueConstant Kt = 0.66 "N.m/A";  
parameter SI.ElectricalTorqueConstant Kb = Kt "N.m/A, equal to Kt for consistent units";  
parameter SI.MomentOfInertia J = 0.0026 "kg.m^2";  
parameter SI.RotationalDampingConstant b = 0.01 "N.m.s/rad";  
parameter SI.Inductance L = 0.05 " H";  

The parameter designation means that the value cannot be changed during simulation. An interesting aspect of Modelica is that the units of these parameters can be specified –SI is a shorthand for “Modelica.SIunits”. Using units is optional, but very useful. If the units had not been specified, all of these parameters would have been of type “Real”. The part in quotes after the value is a comment; C++, Java style of commenting, “//”, can also be used. The units shown in the comment section were written in by the user and, so, must be consistent with the units as defined by Modelica. Having the ability to specify units for physical quantities is part of Modelica’s multi-physical-media modeling capability.

The initial values of the state variables are specified in a similar manner:

SI.Angle theta(start = 0.0) "Shaft angle, rad";  
SI.AngularVelocity omega(start = 0.0) "Shaft angular velocity, rad/s";  
SI.Current i(start = 0.0) "Current in circuit, A";  

The interesting variant here is that for each of the state variables, the initial condition is given as a start value. In this case, the problem is a pure ODE (i.e., a causal problem) so the state variables all have independent initial conditions. If the problem were a DAE (Differential-Algebraic Equation, acausal problem) the initial conditions would not all be independent so the start values would only be suggestions.

Results from simulating this system match the results presented above with Python. As can be seen in this section, for a problem expressed in equation form, there is essentially no “programming” needed; just copy the equations, list the parameter values, and define the system variables.

Modelica — Graphical Modeling Mode

The Modelica Standard Library contains a large set of modeling components from a variety of physical media. This is probably the most common way of creating simulations in Modelica. For the DC motor simulation, this means that components are selected from the library, they are given parameter values, and then the simulation is run. Figure 4 shows the model for the DC motor.

This is a semi-representational model in that some of the components look like their physical counterparts but others can be more abstract. The physical orientation of components on the diagram does not have any significance.

In some sense, the most important component in this diagram is the emf unit. It handles the energy conversion between electrical and rotational-mechanical domains. With the emf unit in place, this model can do proper energy/power bookkeeping across the whole system.

Each of the energy domains is identified by its connectors: blue and white squares for the electrical portion and black and white circles for the rotational section. There is an additional energy domain on the diagram, thermal, identified by red squares, but that is not used in this simulation. A right-click on any of the components will bring up a menu for its parameters. These can be entered directly as numerical values, but it is often preferable to use named parameter variables (as above) and enter the parameter names. That is what is done in this model.

Again, simulation produces the previous results. The modeling at this high level of abstraction can be very efficient if the library contains the needed components. The Modelica Standard Library has a large selection of components so can often be used successfully for simulation of systems crossing multiple energy domains.

Java

As mentioned earlier, Java is a general purpose, object-oriented language and is thus a major abstraction level below either of the Modelica solutions. It requires much more detailed programming (real programming, not equations as in Modelica) in return for which it provides the full power of object-oriented programming and much higher speed.

General purpose languages do not usually include the specialized mathematical tools needed for dynamic system simulation. However, most of these languages have a wide variety of libraries available to fill that gap. For these examples, the Hipparchus library is used to supply the ODE solver 48. As an aside, it is a large library and can meet many other mathematical needs in Java programming. The Hipparchus library is a successor to the Apache Commons Math Library.

As with most ODE solvers, a call-back structure is used where the derivative values of all of the state variables are computed in a function (“method” in Java) and a reference to that function is given to the solver so it can call the derivative function as needed by whatever solution algorithm is being used. In Java, the derivative method is embedded in a class which is what is given to the ODE integration system.

This is what that function looks like for the DC motor system:

public double[] computeDerivatives(double t, double[] y)     
{    
    // ODE: d angle / dt = omega; d omega/dt = (torque - b omega)/J    
    // States: 0 - angle; 1 - omega; 2 - current    
    omega = y[1]; // Angular velocity of shaft    
    ic = y[2]; // Current in circuit    
    vBack = Kb * omega;    
    torq = Kt * ic;    
    double [] dydt = new double[3];    
    // Construct derivatives of states    
    dydt[0] = omega; // d theta / dt = omegga    
    dydt[1] = (torq - b * omega) / J;  // d omega / dt = (torque - damping) / inertia    
    dydt[2] = (v - vBack - R * ic) / L;  // di/dt = (voltage - back-emf - R * current) / L    
    return dydt;     
}    

There are many ways to run a simulation using ODE solvers. In this case, a structure is used that is very convenient for problems where control will be applied. In this structure, the ODE solver is called within an event loop and simulated from one event to the next. This simulation does not actually have any control, although it would be natural to add angular velocity or angular position control. The event that is present is data logging, writing the results to a file for later plotting and/or analysis. The event loop code is:

while(t <= tf)  
{  
    if(t >= (tNextLog - guard))  
    {  
        tNextLog += dtLog;  
        pw.format("%g %g %g %g %g %n", t, y[0], y[1], y[2], volt);  
    }  
    tNextSim = tNextLog; // Simulate until next event  
    ODEStateAndDerivative finalState = dp853.integrate(ode, state, tNextSim);  
    y = finalState.getPrimaryState();  
    t = tNextSim;  
    state = new ODEState(t, y); // setup for next iteration             
}  

Any number of events can be included in a loop like this. The simulation is run until the next event. These events are all time based. Condition-based events (for example, a state variable crossing a boundary value) are trickier and are best handled inside the solver itself, although they can be done iteratively with a structure such as this one. The dp853 refers to the Dormand-Prince-8(5,3) ODE integration algorithm 77. The results are written to a file and plotted using an open-source plotting tool named gnuplot 78. As noted above, one of the advantages of Java is execution speed. By comparison, the Java version is at least 10 times faster than the Modelica-Equations version to compile and execute.