4. Examples

In the next sections, it will be shown how to use the JModelica.org platform for simulation of various FMUs.

The Python commands in these examples may be copied and pasted directly into a Python shell, in some cases with minor modifications. Alternatively, they may be copied into a text file, which also is the recommended way.

4.1. Simulation of a high-index model

Mechanical component-based models often result in high-index DAEs. In order to efficiently integrate such models, Modelica tools typically employs an index reduction scheme, where some equations are differentiated, and dummy derivatives are selected. In order to demonstrate this feature, we consider the model Modelica.Mechanics.Rotational.Examples.First from the Modelica Standard library, see Figure 5.2, “Modelica.Mechanics.Rotational.First connection diagram”. The model is of high index since there are two rotating inertias connected with a rigid gear.

Figure 5.2. Modelica.Mechanics.Rotational.First connection diagram

Modelica.Mechanics.Rotational.First connection diagram

First create a Python script file and enter the usual imports:

import matplotlib.pyplot as plt
from pymodelica import compile_fmu
from pyfmi import load_fmu

Next, the model is compiled and loaded:

# Compile model
fmu_name = compile_fmu("Modelica.Mechanics.Rotational.Examples.First",())

# Load model
model = load_fmu(fmu_name)

Notice that no file name, just an empty tuple, is provided to the function compile_fmu, since in this case the model that is compiled resides in the Modelica standard library. In the compilation process, the index reduction algorithm is invoked. Next, the model is simulated for 3 seconds:

# Load result file
res = model.simulate(final_time=3.)

Finally, the simulation results are retrieved and plotted:

w1 = res['inertia1.w']
w2 = res['inertia2.w']
w3 = res['inertia3.w']
tau = res['torque.tau']
t = res['time']

plt.figure(1)
plt.subplot(2,1,1)
plt.plot(t,w1,t,w2,t,w3)
plt.grid(True)
plt.legend(['inertia1.w','inertia2.w','inertia3.w'])
plt.subplot(2,1,2)
plt.plot(t,tau)
plt.grid(True)
plt.legend(['tau'])
plt.xlabel('time [s]')
plt.show()

You should now see a plot as shown below.

Figure 5.3. Simulation result for Modelica.Mechanics.Rotational.Examples.First

Simulation result for Modelica.Mechanics.Rotational.Examples.First

4.2. Simulation and parameter sweeps

This example demonstrates how to run multiple simulations with different parameter values. Sweeping parameters is a useful technique for analysing model sensitivity with respect to uncertainty in physical parameters or initial conditions. Consider the following model of the Van der Pol oscillator:

  model VDP
    // State start values
    parameter Real x1_0 = 0;
    parameter Real x2_0 = 1;

    // The states
    Real x1(start = x1_0);
    Real x2(start = x2_0);

    // The control signal
    input Real u;

  equation
    der(x1) = (1 - x2^2) * x1 - x2 + u;
    der(x2) = x1;
  end VDP;

Notice that the initial values of the states are parametrized by the parameters x1_0 and x2_0. Next, copy the Modelica code above into a file VDP.mo and save it in your working directory. Also, create a Python script file and name it vdp_pp.py. Start by copying the commands:

import numpy as N
import pylab as P
from pymodelica import compile_fmu
from pyfmi import load_fmu

into the Python file. Compile and load the model:

# Define model file name and class name
model_name = 'VDP'
mofile = 'VDP.mo'

# Compile model
fmu_name = compile_fmu(model_name,mofile)

Next, we define the initial conditions for which the parameter sweep will be done. The state x2 starts at 0, whereas the initial condition for x1 is swept between -3 and 3:

# Define initial conditions
N_points = 11
x1_0 = N.linspace(-3.,3.,N_points)
x2_0 = N.zeros(N_points)

In order to visualize the results of the simulations, we open a plot window:

fig = P.figure()
P.clf()
P.hold(True)
P.xlabel('x1')
P.ylabel('x2')

The actual parameter sweep is done by looping over the initial condition vectors and in each iteration set the parameter values into the model, simulate and plot:

for i in range(N_points):
    # Load model
    vdp = load_fmu(fmu_name)  
    # Set initial conditions in model
    vdp.set('x1_0',x1_0[i])
    vdp.set('x2_0',x2_0[i])
    # Simulate 
    res = vdp.simulate(final_time=20)
    # Get simulation result
    x1=res['x1']
    x2=res['x2']
    # Plot simulation result in phase plane plot
    P.plot(x1, x2,'b')
P.grid()
P.show()

You should now see a plot similar to that in Figure 5.4, “Simulation result-phase plane”.

Figure 5.4. Simulation result-phase plane

Simulation result-phase plane

4.3. Simulation of an Engine model with inputs

In this example the model is larger than the previous. It is a slightly modified version of the model EngineV6_analytic from the Multibody library in the Modelica Standard Library. The modification consists of a replaced load with a user defined load. This has been done in order to be able to demonstrate how inputs are set from a Python script. In Figure 5.5, “Overview of the Engine model” the model is shown.

Figure 5.5. Overview of the Engine model

Overview of the Engine model

The Modelica code for the model is shown below, copy and save the code in a file named EngineV6.mo.

model EngineV6_analytic_with_input
  output Real engineSpeed_rpm= Modelica.SIunits.Conversions.to_rpm(load.w);
  output Real engineTorque = filter.u;
  output Real filteredEngineTorque = filter.y;
  
  input Real u;
  
  import Modelica.Mechanics.*;

  inner MultiBody.World world;
  MultiBody.Examples.Loops.Utilities.EngineV6_analytic engine(redeclare 
      model Cylinder = MultiBody.Examples.Loops.Utilities.Cylinder_analytic_CAD);
      
  Rotational.Components.Inertia load(
    phi(start=0,fixed=true), w(start=10,fixed=true),
    stateSelect=StateSelect.always,J=1);
  Rotational.Sensors.TorqueSensor torqueSensor;
  Rotational.Sources.Torque torque;
  
  Modelica.Blocks.Continuous.CriticalDamping filter(
    n=2,initType=Modelica.Blocks.Types.Init.SteadyState,f=5);

equation 
  torque.tau = u;
  
  connect(world.frame_b, engine.frame_a);
  connect(torque.flange, load.flange_b);
  connect(torqueSensor.flange_a, engine.flange_b);
  connect(torqueSensor.flange_b, load.flange_a);
  connect(torqueSensor.tau, filter.u);
  annotation (experiment(StopTime=1.01));

end EngineV6_analytic_with_input;

Now that the model has been defined, we create our Python script which will compile, simulate and visualize the result for us. Create a new text-file and start by copying the below commands into the file. The code will import the necessary methods and packages into Python.

from pymodelica import compile_fmu
from pyfmi import load_fmu
import pylab as P

Compiling the model is performed by invoking the compile_fmu method where the first argument is the name of the model and the second argument is where the model is located (which file). The method will create an FMU in the current directory and in order to simulate the FMU, we need to additionally load the created FMU into Python. This is done with the load_fmu method which takes the name of the FMU as input.

name = compile_fmu("EngineV6_analytic_with_input", "EngineV6.mo")

model = load_fmu(name)

So, now that we have compiled the model and loaded it into Python we are almost ready to simulate the model. First however, we retrieve the simulation options and specify how many result points we want to receive after a simulation.

opts = model.simulate_options()
opts["ncp"] = 1000 #Specify that 1000 output points should be returned

A simulation is finally performed using the simulate method on the model and as we have changed the options, we need to additionally provide these options to the simulate method.

res = model.simulate(options=opts)

The simulation result is returned and stored into the res object. Result for a trajectory is easily retrieved using a Python dictionary syntax. Below is the visualization code for viewing the engine torque and the engine speed. One could instead use the Plot GUI for the visualization as the result are stored in a file in the current directory.

P.subplot(211)
P.suptitle("EngineV6")
P.plot(res["time"],res["filteredEngineTorque"], label="Filtered Engine Torque")
P.grid()
P.legend()
P.ylabel("Torque [N.m]")
P.subplot(212)
P.plot(res["time"],res["engineSpeed_rpm"], label="Engine Speed")
P.grid()
P.legend()
P.xlabel("Time [s]")
P.ylabel("Speed [1/min]")
P.show()

In Figure 5.6, “Resulting trajectories for the engine model.” the trajectories are shown.

Figure 5.6. Resulting trajectories for the engine model.

Resulting trajectories for the engine model.

Above we have simulated the engine model and looked at the result, we have not however specified any load as input. Remember that the model we are looking at has a user specified load. Now we will create a Python function that will act as our input. We create a function that depends on the time and returns the value for use as input.

def input_func(t):
    return -100.0*t

In order to use this input in the simulation, simply provide the name of the input variable and the function as the input argument to the simulate method, see below.

res = model.simulate(options=opts, input=("u",input_func))

Simulate the model again and look at the result and the impact of the input.

Large models contain an enormous amount of variables and by default, all of these variables are stored in the result. Storing the result takes time and for large models the saving of the result may be responsible for the majority of the overall simulation time. Not all variables may be of interest, for example in our case, we are only interested in two variables so storing the other variables are not necessary. In the options dictionary there is a filter option which allows to specify which variables should be stored, so in our case, try the below filter and look at the impact on the simulation time.

opts["filter"] = ["filteredEngineTorque", "engineSpeed_rpm"]

4.4. Simulation using the native FMI interface

This example shows how to use the native JModelica.org FMI interface for simulation of an FMU (FMI 1.0 for Model Exchange). The FMU that is to be simulated is the bouncing ball example from Qtronics FMU SDK (http://www.qtronic.de/en/fmusdk.html). This example is written similar to the example in the documentation of the 'Functional Mock-up Interface for Model Exchange' version 1.0 (https://www.fmi-standard.org/). The bouncing ball model is to be simulated using the explicit Euler method with event detection.

The example can also be found in the Python examples catalog in the JModelica.org platform.

The bouncing ball consists of two equations,

and one event function (also commonly called root function),

Where the ball bounces and lose some of its energy according to,

Here, h is the height, g the gravity, v the velocity and e a dimensionless parameter. The starting values are, h=1 and v=0 and for the parameters, e=0.7 and g = 9.81.

4.4.1. Implementation

Start by importing the necessary modules,

import numpy as N 
import pylab as P                  # Used for plotting
from pyfmi.fmi import FMUModelME1  # The FMI Interface for Model Exchange

Next, the FMU is to be loaded and initialized

# Load the FMU by specifying the fmu together with the path.
bouncing_fmu = FMUModelME1('/path/to/FMU/bouncingBall.fmu')

Tstart = 0.5                 # The start time.
Tend   = 3.0                 # The final simulation time.  
bouncing_fmu.time = Tstart   # Set the start time before the initialization.
                             # (Defaults to 0.0)
bouncing_fmu.initialize()    # Initialize the model. Also sets all the start 
                             # attributes defined in the XML file.

The first line loads the FMU and connects the C-functions of the model to Python together with loading the information from the XML-file. The start time also needs to be specified by setting the property time. The model is also initialized, which must be done before the simulation is started.

Note that if the start time is not specified, FMUModelME1 tries to find the starting time in the XML-file structure 'default experiment' and if successful starts the simulation from that time. Also if the XML-file does not contain any information about the default experiment the simulation is started from time zero.

Then information about the first step is retrieved and stored for later use.

# Get Continuous States
x = bouncing_fmu.continuous_states
# Get the Nominal Values
x_nominal = bouncing_fmu.nominal_continuous_states
# Get the Event Indicators
event_ind = bouncing_fmu.get_event_indicators()
        
# Values for the solution
vref  = [bouncing_fmu.get_variable_valueref('h')] + \
        [bouncing_fmu.get_variable_valueref('v')]        # Retrieve the valureferences for the
                                                         # values 'h' and 'v'
t_sol = [Tstart]
sol = [bouncing_fmu.get_real(vref)]

Here the continuous states together with the nominal values and the event indicators are stored to be used in the integration loop. In our case the nominal values are all equal to one. This information is available in the XML-file. We also create lists which are used for storing the result. The final step before the integration is started is to define the step-size.

time = Tstart
Tnext = Tend   # Used for time events
dt = 0.01      # Step-size

We are now ready to create our main integration loop where the solution is advanced using the explicit Euler method.

# Main integration loop.
while time < Tend and not bouncing_fmu.get_event_info().terminateSimulation:
    #Compute the derivative of the previous step f(x(n), t(n))
    dx = bouncing_fmu.get_derivatives()
    
    # Advance
    h = min(dt, Tnext-time)
    time = time + h
            
    # Set the time
    bouncing_fmu.time = time
            
    # Set the inputs at the current time (if any)
    # bouncing_fmu.set_real,set_integer,set_boolean,set_string (valueref, values)
            
    # Set the states at t = time (Perform the step using x(n+1)=x(n)+hf(x(n), t(n))
    x = x + h*dx 
    bouncing_fmu.continuous_states = x

This is the integration loop for advancing the solution one step. The loop continues until the final time have been reached or if the FMU reported that the simulation is to be terminated. At the start of the loop the derivatives of the continuous states are retrieved and then the simulation time is incremented by the step-size and set to the model. It could also be the case that the model depends on inputs which can be set using the set_(real/...) methods.

Note that only variables defined in the XML-file to be inputs can be set using the set_(real/...) methods according to the FMI specification.

The step is performed by calculating the new states (x+h*dx) and setting the values into the model. As our model, the bouncing ball also consist of event functions which needs to be monitored during the simulation, we have to check the indicators which is done below.

    # Get the event indicators at t = time
    event_ind_new = bouncing_fmu.get_event_indicators()
            
    # Inform the model about an accepted step and check for step events
    step_event = bouncing_fmu.completed_integrator_step()
            
    # Check for time and state events
    time_event  = abs(time-Tnext) <= 1.e-10
    state_event = True if True in ((event_ind_new>0.0) != (event_ind>0.0)) else False

Events can be, time, state or step events. The time events are checked by continuously monitoring the current time and the next time event (Tnext). State events are checked against sign changes of the event functions. Step events are monitored in the FMU, in the method completed_integrator_step() and return True if any event handling is necessary. If an event have occurred, it needs to be handled, see below.

    # Event handling
    if step_event or time_event or state_event:
        eInfo = bouncing_fmu.get_event_info()
        eInfo.iterationConverged = False
                
        # Event iteration
        while eInfo.iterationConverged == False:
            bouncing_fmu.event_update('0')       # Stops at each event iteration
            eInfo = bouncing_fmu.get_event_info()

            # Retrieve solutions (if needed)
            if eInfo.iterationConverged == False:
                # bouncing_fmu.get_real,get_integer,get_boolean,get_string(valueref)
                pass
                
        # Check if the event affected the state values and if so sets them
        if eInfo.stateValuesChanged:
            x = bouncing_fmu.continuous_states
            
        # Get new nominal values.
        if eInfo.stateValueReferencesChanged:
            atol = 0.01*rtol*bouncing_fmu.nominal_continuous_states
                    
        # Check for new time event
        if eInfo.upcomingTimeEvent:
            Tnext = min(eInfo.nextEventTime, Tend)
        else:
            Tnext = Tend

If an event occurred, we enter the iteration loop where we loop until the solution of the new states have converged. During this iteration we can also retrieve the intermediate values with the normal get methods. At this point eInfo contains information about the changes made in the iteration. If the state values have changed, they are retrieved. If the state references have changed, meaning that the state variables no longer have the same meaning as before by pointing to another set of continuous variables in the model, for example in the case with dynamic state selection, new absolute tolerances are calculated with the new nominal values. Finally the model is checked for a new time event.

    event_ind = event_ind_new
        
    # Retrieve solutions at t=time for outputs
    # bouncing_fmu.get_real,get_integer,get_boolean,get_string (valueref)
            
    t_sol += [time]
    sol += [bouncing_fmu.get_real(vref)]

In the end of the loop, the solution is stored and the old event indicators are stored for use in the next loop.

After the loop have finished, by reaching the final time, we plot the simulation results

# Plot the height
P.figure(1)
P.plot(t_sol,N.array(sol)[:,0])
P.title(bouncing_fmu.get_name())
P.ylabel('Height (m)')
P.xlabel('Time (s)')
# Plot the velocity
P.figure(2)
P.plot(t_sol,N.array(sol)[:,1])
P.title(bouncing_fmu.get_name())
P.ylabel('Velocity (m/s)')
P.xlabel('Time (s)')
P.show()

and the figure below shows the results.

Figure 5.7. Simulation result

Simulation result

4.5. Simulation of Co-Simulation FMUs

Simulation of a Co-Simulation FMU follows the same workflow as simulation of a Model Exchange FMU. The model we would like to simulate is a model of a bouncing ball, the file bouncingBall.fmu is located in the examples folder in the JModelica.org installation, pyfmi/examples/files/CS1.0/. The FMU is a Co-simulation FMU and in order to simulate it, we start by importing the necessary methods and packages into Python:

import pylab as P           # For plotting
from pyfmi import load_fmu  # For loading the FMU 

Here, we have imported packages for plotting and the method load_fmu which takes as input an FMU and then determines the type and returns the appropriate class. Now, we need to load the FMU.

model = load_fmu('bouncingBall.fmu')

The model object can now be used to interact with the FMU, setting and getting values for instance. A simulation is performed by invoking the simulate method:

res = model.simulate(final_time=2.)

As a Co-Simulation FMU contains its own integrator, the method simulate calls this integrator. Finally, plotting the result is done as before:

# Retrieve the result for the variables
h_res = res['h']
v_res = res['v']
t     = res['time'] 
# Plot the solution
# Plot the height
fig = P.figure()
P.clf()
P.subplot(2,1,1)
P.plot(t, h_res)
P.ylabel('Height (m)')
P.xlabel('Time (s)')
# Plot the velocity
P.subplot(2,1,2)
P.plot(t, v_res)
P.ylabel('Velocity (m/s)')
P.xlabel('Time (s)')
P.suptitle('FMI Bouncing Ball')
P.show()

and the figure below shows the results.

Figure 5.8. Simulation result

Simulation result