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.
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.
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.
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
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”.
In this example the model is rather larger then 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.
The Modelica code for the model is shown below, copy and save the code in a file named
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.
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"]
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.
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 is depended on inputs which can be set using the
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 monitor 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.
Simulation of an Co-Simulation FMU follow the same workflow as simulation of an 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')
model object can now be used to interact with the FMU, setting and getting values for instance. A simulation is performed by invoking the
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.