Volume 10, Issue 1

Articles
Trott's Corner
New Products
New Publications
Calendar
News Bulletins
New Resources
Classifieds

Editorial Policy
Staff and Contributors
Submissions
Subscriptions
Back Issues
Contact Information

MathModelica--An Object-Oriented Mathematical Modeling and Simulation Environment

# 4. Application Examples

This section gives a number of application examples of the use of the MathModelica environment. The intent is to demonstrate the power of integration and interactivity--the interplay between the object-oriented modeling and simulation capabilities of Modelica integrated with the powerful scripting facilities of Mathematica within MathModelica. This includes the representation of simulation results as 1D and 2D interpolating functions of time being combined with arithmetic operations and functions in expressions, advanced plotting facilities, and computational capabilities, such as design optimization, fourier analysis, and solution of time-dependent PDEs.

## 4.1. Advanced Plotting and Interpolating Functions

This section illustrates the flexible usage of simulation results (represented as interpolating functions) for further computations that may include simulation results in expressions, and for both simple and advanced plotting. The following simple bouncing ball model [12] is used in the simulation and plotting examples.

### Interpolating Function Representation of Simulation Results

The following simulation of the previous BouncingBall model is done for a short time period using very few points.

The results returned by Simulate are represented by an access descriptor or handle. Note that the output also mentions the parameters c and g in the "variables" list even though their values are constant and not generated by the simulation. Some of the contents of such a descriptor are shown as the result of the previous call to Simulate. At this stage the simulation data is stored on disk and referenced by res1, which acts as a handle to the simulation data. When one of the variables from the last simulation is referenced, e.g. height, radius, etc., the data for that variable are loaded into the system in a load-by-need manner and represented as an ordinary Mathematica InterpolatingFunction.

Working with simulation result datasets in Mathematica is done in a very convenient way using the Mathematica InterpolatingFunction mechanism that encapsulates the data into a function object and provides interpolation so that the data acts as a regular function. The mechanism of loading simulation data into the system and representing it as a function object is performed by the function VariableTrajectory, which can be called explicitly as follows, but is called automatically on any variable from the last simulation when referenced.

The expression returned from the VariableTrajectory is an anonymous function object having one formal parameter (for the time t) and consists of a body containing an expression that computes a value of the variable for the given time. In this case the body consists of a Mathematica Which expression that switches between two (or more) InterpolatingFunction objects, dependent on whether the time is less than the event time point at 0.428 or not. The InterpolatingFunction objects uses interpolation order 3 as the default but can be altered by using options. Pure discrete data (i.e, data changing only at event points) is encapsulated by one InterpolatingFunction object with zero interpolation order to get a piecewise constant behavior in the interpolation. This is more efficient than using Which statements. The system will automatically choose the most efficient representation of these two alternatives.

The interpolation function can now be used in any computation in Mathematica. In this case we just evaluate the derivative at the time 0.2.

In the previous case Mathematica made the differentiation of the InterpolatingFunction object h. Normally the derivatives of the simulation variables are also available in the simulation data.

As previously mentioned, to help the MathModelica user, variables of the most recent simulation are always accessible directly. In this case the function VariableTrajectory is automatically applied. Therefore, instead of assigning the variable h as done earlier, we can write the following and get the same result.

Note, to keep variable values from previous simulations accessible, we should use VariableTrajectory on the appropriate variable, specifying the desired descriptor (e.g., res1).

Now we perform a new simulation, with the result denoted by res2.

Having the previous height curve represented as the function object h, we can easily compute the difference of the curves between the simulations (e.g., using a plot expression height[t]-h[t]).

### PlotSimulation

First, we simulate the bouncing ball for eight seconds and store the results in the variable res1 for subsequent use in the plotting examples.

The command PlotSimulation is used for simple standard plots. If nothing else is specified (i.e., by the optional SimulationResult parameter), the command refers to the results from the most recent simulation. In the following diagram the height above ground of the ball from the bouncing ball model simulation is plotted for the first eight seconds of simulation. The optional parameter PlotJoined has been set to False to create a dotted plot.

Figure 16. Dotted plot of bouncing ball example model.

PlotSimulation can also handle expressions containing simulated results. When this is done, a warning is returned to emphasize that interpolation is performed, which could result in a slightly less accurate plot.

Figure 17. Plot of expression involving interpolated function of simulation result.

Plotting several arbitrary functions can be done using a list of function expressions instead of a single expression.

Figure 18. Plot of arbitrary functions in the same diagram.

Now we simulate the bouncing ball again but with a different value of the coefficient of restitution, c, which is changed to 0.95. The result is stored in res2.

The optional argument SimulationResult specifies which simulation data to use. In this case we will use res1 and res2. The two plots are stored as two graphics objects gr1 and gr2, which are displayed together using the Mathematica command Show.

Figure 19. Parallel display of two diagrams.

It is possible to plot variables with the same name from several different simulations together. This is specified by an array value for the optional argument SimulationResult.

Figure 20. Plot of variables from several simulations in the same diagram.

The plot colors are specified by the \$PlotSimulationColors predefined variable.

Here we check the color scheme of \$PlotSimulationColors.

Figure 21. Plot of multiple curves with different colors.

### ParametricPlotSimulation

Parametric plots can be done using ParametricPlotSimulation.

Figure 22. A parametric plot.

In the same way as PlotSimulation, the ParametricPlotSimulation function can handle several plots.

Figure 23. Multiple parametric plots in the same diagram.

ParametricPlotSimulation can also handle results from different simulations and plot only the actual data points.

Figure 24. Parametric plots of data points from different simulations in the same diagram.

### ParametricPlotSimulation3D

In this example we are going to use the Rossler attractor to show the ParametricPlotSimulation3D command. The Rossler attractor is named after Otto Rossler's work in chemical kinetics. The system is described by three coupled nonlinear differential equations.

Here and are constants. The attractor never forms limit circles nor does it ever reach a steady state. The model is shown in Mathematica syntax, enabling the use of Greek characters.

The model is simulated using different initial values. Changing these can considerably influence the appearance of the attractor.

The Rossler attractor is easy to plot using ParametricPlotSimulation3D:

Figure 25. 3D parametric plot of interpolated curve from the Rossler attractor simulation.

The plot does not look smooth at some areas, especially for high values of Z. Let us take a look at the actual data points of the simulation.

Figure 26. 3D parametric plot of actual data points from the Rossler attractor simulation.

There seem to be few data points outside the "rings." This can be fixed by adding data points during simulation. The default value is 500. Let us try 1000 data points.

Now the plot looks smoother.

Figure 27. 3D parametric plot of a curve with many data points from the Rossler attractor simulation.

## 4.2. Design Optimization

Let us describe an example of how the powerful Mathematica scripting language available within MathModelica can be utilized to solve nontrivial optimization problems that contain dynamic simulations. First, we will define a Modelica model of a linear actuator with spring damped stopping and then a first-order system. Using MathModelica scripting, we will then find a damping for the translational spring-damper such that the step response is as "close" as possible to the step response from a first-order system.

Consider the following model of a linear actuator with a spring damped connection to an anchoring point.

Figure 28. A LinearActuator model containing a spring damped connection to an anchoring point.

Here is the corresponding model code in Mathematica-style syntax (the squares contain graphical annotations for icons, lines, and so on).

As a comparison, we also show the model code in Modelica textual syntax:

The small square boxes in the code represent graphical annotations about the visual appearance of the components of the model, as seen in Figure 28. These annotations are hidden behind boxes; otherwise, the code becomes harder to read.

Here we simulate a step response and store the result in res0.

Figure 29. Plot of the step response from the linear actuator.

Assume that we have some freedom in choosing the damping in the translational spring-damper. A number of simulation runs shows what kind of behavior we have for different values of the damping parameter d. The Mathematica Table function is used with Simulate to collect the results into an array res. This array then contains the results from simulations of LinearActuator with a damping of 2 to 14 with a step size of 2 (i.e., seven simulations are performed).

Figure 30. Plots of step responses from seven simulations of the linear actuator with different camping coefficients.

Now assume that we would like to choose the damping d so that the resulting system behaves as closely as possible to the following first-order system response, obtained by solving a first-order ODE using NDSolve:

Clear[y]

We make a comparison with the step response we simulated first (d=2) and the first-order system.

Figure 31. Comparison plot between the step response of the linear actuator and a first-order system.

Now, let us make things a little more automatic. Simulate and compute the integral of the square error from to .

Here we define a function, f[a_], doing the same thing as previously, but for different spring-damper parameter values

and tabulate some results for .

The tabulated values are interpolated using an interpolating function object. The default interpolation order is 3.

Figure 32. Plot of the error function for finding a minimum deviation from the desired step response.

The minimizing value of can be computed using FindMinimum.

Here is a simulation with the optimal parameter value.

Here is a plot comparing the first- and second-order system response with a plot of the squared error amplified by a factor of 100.

Figure 33. Comparison plot of the first- and second-order system step responses with the squared error.

## 4.3. Fourier Analysis of Simulation Data

Consider a weak axis excited by a torque pulse train. The axis is modeled by three segments joined by two torsion springs. The following diagram is imported from the MathModelica Model Editor where the model was defined.

Figure 34. A WeakAxis model excited by a torque pulse train.

Here is the corresponding Modelica code in Mathematica syntax.

The model code in Modelica syntax.

Here we simulate the model during 200 seconds.

The plot of the angular velocity of the rightmost axis segment appears as follows.

Figure 35. Plot of the angular velocity of the rightmost axis segment of the WeakAxis model.

Now, let us sample the interpolated function Inertia3.w using a sample frequency of 4Hz, and put the result into an array using the Mathematica Table array constructor.

Here we compute the absolute values of the discrete Fourier transform of data1 with the mean value removed.

Here we plot the first 80 points of the data.

Figure 36. Plot of the data points of the Fourier transformed angular velocity.

It is easy to write a function FourierPlot that repeats the previous operations. FourierPlot also scales the axes such that amplitudes of trigonometric components are plotted against frequency (Hz).

Figure 37. Plot of the curve of the Fourier transformed angular velocity.

It can be shown that the frequencies of the eigenmodes of the system are given by the imaginary parts of the eigenvalues of the following matrix ( and are the spring constants).

These values, 0.256077, 0.143344, fit very well with the peaks in the previous diagram.

## 4.4. Solution and 2D-Interpolation of Discretized PDEs

Currently Modelica cannot handle partial differential equations directly since there is only the notion of differentiation with respect to time built into the language. However, in many cases derivatives with respect to other variables such as, for example, spatial dimensions can be handled by simple discretizations schemes easily implemented using the array capabilities in Modelica. Here we will give an example of how the one-dimensional wave equation can be represented in Modelica and how MathModelica can be used for simulation and display of the results, as well as representing the result as a two-dimensional interpolating function.

The one-dimensional wave equation is given by a partial differential equation of the following form

where is a function of both space and time. As a physical example, let us consider a duct of length 10 where we let describe its spatial dimension. Since Modelica can only handle time as the independent variable, we need to discretize the problem in the spatial dimension and approximate the spatial derivatives using difference approximations using the approximation

Utilizing this approach a MathModelica model of a duct whose pressure dynamics is given by the wave equation can be written as follows.

Here we are using a Modelica function initialPressure (defined as follows) to specify the initial value of the pressure along the duct. Assume that we would like an initial pressure profile in the duct of the form (Figure 38).

Figure 38. The initial pressure profile of the duct.

To provide an initial value profile (the start value of the vector variable p) for the pressure, we need to define a function that returns a sampled version of this profile. In principle this function could have been expressed in ordinary Mathematica function syntax with type extensions as in the MathCode system [41], but here we are using a more Modelica-like Mathematica syntax that is close to the syntax used for declaration of classes.

Here we simulate the WaveEquationSample model.

The result is packed into a 2D interpolation function.

A 3D plot of the pressure distribution in the duct can be obtained as follows.

Figure 39. A 3D plot of the pressure distribution in the duct.

Let us also plot the wave as a function of position for a fixed point in time.

Figure 40. The wave as a function of position at time=1.2.