This example will demonstrate how a model with two inputs with data from a matlab-file can be simulated. The model to be simulated is a quadruple tank connected to two pumps, which is also the inputs to the model. The model is depicted in Figure 7.2 and in the code below the corresponding Modelica code is listed.

model QuadTank
    // Process parameters
  parameter Modelica.SIunits.Area A1=4.9e-4, A2=4.9e-4, A3=4.9e-4, A4=4.9e-4;
  parameter Modelica.SIunits.Area a1=0.03e-4, a2=0.03e-4, a3=0.03e-4, a4=0.03e-4;
  parameter Modelica.SIunits.Acceleration g=9.81;
  parameter Real k1_nmp(unit="m3^/s/V") = 0.56e-6, k2_nmp(unit="m^3/s/V") = 0.56e-6;
  parameter Real g1_nmp=0.30, g2_nmp=0.30;

   // Initial tank levels
  parameter Modelica.SIunits.Length x1_0 = 0.06270;
  parameter Modelica.SIunits.Length x2_0 = 0.06044;
  parameter Modelica.SIunits.Length x3_0 = 0.02400;
  parameter Modelica.SIunits.Length x4_0 = 0.02300;

   // Tank levels
  Modelica.SIunits.Length x1(start=x1_0,min=0.0001/*,max=0.20*/);
  Modelica.SIunits.Length x2(start=x2_0,min=0.0001/*,max=0.20*/);
  Modelica.SIunits.Length x3(start=x3_0,min=0.0001/*,max=0.20*/);
  Modelica.SIunits.Length x4(start=x4_0,min=0.0001/*,max=0.20*/);

  // Inputs
  input Modelica.SIunits.Voltage u1;
  input Modelica.SIunits.Voltage u2;

 equation
   der(x1) = -a1/A1*sqrt(2*g*x1) + a3/A1*sqrt(2*g*x3) +
                                   g1_nmp*k1_nmp/A1*u1;
   der(x2) = -a2/A2*sqrt(2*g*x2) + a4/A2*sqrt(2*g*x4) +
                                   g2_nmp*k2_nmp/A2*u2;
   der(x3) = -a3/A3*sqrt(2*g*x3) + (1-g2_nmp)*k2_nmp/A3*u2;
   der(x4) = -a4/A4*sqrt(2*g*x4) + (1-g1_nmp)*k1_nmp/A4*u1;

end QuadTank;

Lets begin with the the example, copy and paste the Modelica code and save it into QuadTank.mo and open a python script file. We start by importing the necessary objects:

from scipy.io.matlab.mio import loadmat
import matplotlib.pyplot as plt
import numpy as N

from pymodelica import compile_jmu
from pyjmi import JMUModel

The input data is stored in qt_par_est_data.mat which can be found in the examples/files catalogue in JModelica.org. Copy it into your working directory and paste the following commands to load the data-file and extract the data trajectories:

data = loadmat('qt_par_est_data.mat',appendmat=False)

# Extract data series  
t_meas = data['t'][6000::100,0]-60  
u1 = data['u1_d'][6000::100,0]
u2 = data['u2_d'][6000::100,0]

The trajectories have now been extracted and needs to be stacked into a data matrix with the first column as the time vector and the following columns the input of u1 and u2. The names of the variables needs also be connected in the input object:

# Build input trajectory matrix for use in simulation
u_data = N.transpose(N.vstack((t_meas,u1,u2)))
input_object = (['u1','u2'], u_data)

Next, we compile and load the model:

# compile JMU
jmu_name = compile_jmu('QuadTank', 'QuadTank.mo')

# Load model
model = JMUModel(jmu_name)

Now, that the model is compiled and the input have been adapted, lets give the information to the simulate method and simulate:

# Simulate model with input trajectories
res = model.simulate(final_time=60, input=input_object)

The result is retrieved by accessing the res variable as a dictionary with the variable name as key:

x1_sim = res['x1']
x2_sim = res['x2']
x3_sim = res['x3']
x4_sim = res['x4']
u1_sim = res['u1']
u2_sim = res['u2']
t_sim  = res['time']

And then plotted with the help from matplotlib:

plt.figure(1)
plt.subplot(2,2,1)
plt.plot(t_sim,x3_sim)
plt.title('x3')
plt.subplot(2,2,2)
plt.plot(t_sim,x4_sim)
plt.title('x4')
plt.subplot(2,2,3)
plt.plot(t_sim,x1_sim)
plt.title('x1')
plt.xlabel('t[s]')
plt.subplot(2,2,4)
plt.plot(t_sim,x2_sim)
plt.title('x2')
plt.xlabel('t[s]')
plt.show()

plt.figure(2)
plt.subplot(2,1,1)
plt.plot(t_sim,u1_sim,'r')
plt.title('u1')
plt.subplot(2,1,2)
plt.plot(t_sim,u2_sim,'r')
plt.title('u2')
plt.xlabel('t[s]')
plt.show()

In Figure 7.3 the result of the tank levels are shown and in Figure 7.4 the input signals are shown.

The model which is to be simulated is an electric circuit. The model is depicted in Figure 7.5 and consists of resistances, inductors and a capacitor. The circuit is connected to a voltage source which generates a square-wave with an amplitude of 1.0 and a frequency of 0.6 Hz. The model is also available from the examples in the file RLC_Circuit.mo.

This examples assumes that the file RLC_Circuit.mo is located in the working directory.

Start by creating a Python script file and write or (copy paste) the command for importing the model object and for compiling a model together with the library used for plotting:

# Import the function for compilation of models and the JMUModel class
from pymodelica import compile_jmu
from pyjmi import JMUModel

# Import the plotting library
import matplotlib.pyplot as plt

Next, we compile and load the model:

# Compile model
jmu_name = compile_jmu("RLC_Circuit_Square","RLC_Circuit.mo")

# Load model
rlc = JMUModel(jmu_name)

Now we are ready to simulate our model. We are interested in simulating the model from 0.0 to 20.0 seconds. The start time is default to 0.0 so no need to change that, but the final time needs to be changed:

res = rlc.simulate(final_time=20.0) #Simulate the model from 0.0 to 20.0 seconds

After a successful simulation the statistics are printed in the prompt and the results are stored in the variable 'res'. To view the result, we have to retrieve information about the variables we are interested of which is easily done in the following way:

square_y    = res['square.y']
resistor_v  = res['resistor.v']
inductor1_i = res['inductor1.i']
time        = res['time']

And then plotted with the help from matplotlib,

plt.figure(1)
plt.plot(time, square_y, time, resistor_v, time, inductor1_i)    
plt.legend(('square.y','resistor.v','inductor1.i'))
plt.show()

The simulation result is shown in Figure 7.6.

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 7.7. 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 FMUModel

Next, the model is compiled and loaded:

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

# Load model
model = FMUModel(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 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 FMUModel

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)

# Load model
vdp = FMUModel(fmu_name)

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):
    # 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 7.9.

This example will show how to use JModelica.org to simulate an Optimica model and calculate sensitivities of the state variables with respect to a number of free parameters.

The model equations is taken from the Robertson example in the Sundials suite (https://computation.llnl.gov/casc/sundials/main.html) and the model is shown in the code below.

optimization Robertson
    parameter Real p1(free=true)=0.040;
    parameter Real p2(free=true)=1.0e4;
    parameter Real p3(free=true)=3.0e7;
    
    Real y1(start=1.0, fixed=true);
    Real y2(start=0.0, fixed=true);
    Real y3(start=0.0);
  equation
    der(y1) = -p1*y1 + p2*y2*y3;
    der(y2) =  p1*y1 - p2*y2*y3 - p3*(y2*y2);
    0.0 = y1 + y2 + y3 - 1;
end Robertson;

In the model we have set the parameters to free which means that we want to calculate sensitivities of the states with respect to the free parameters.

Lets begin with the the example, copy and paste the Optimica code and save it into Robertson.mop and open a python script file. We start by importing the necessary objects:

# Import the function for compilation of models and the JMUModel class
from pymodelica import compile_jmu
from pyjmi import JMUModel

# Import the plotting library
import matplotlib.pyplot as plt

Next, we compile and load the model:

# Compile model
jmu_name = compile_jmu("Robertson","Robertson.mop")

# Load model
model = JMUModel(jmu_name)

Notice that sensitivity computations are only supported for JMUModels currently. Now that the model is loaded, we have to change the option to activate the sensitivity calculations and also the absolute tolerances:

# Get and set the options
opts = model.simulate_options() #Get the options
opts['IDA_options']['atol'] = [1.0e-8, 1.0e-14, 1.0e-6] #Change the tolerance
opts['IDA_options']['sensitivity'] = True #Activate the sensitivity calculations
opts['ncp'] = 400 #Change the number of communication points

Now lets simulate the model:

res = model.simulate(final_time=4, options=opts)

The sensitivity results are stored as d{variable name}/d{parameter name} in the result object. We are interested in the following sensitivities:

dy1dp1 = res['dy1/dp1']
dy2dp1 = res['dy2/dp1']
dy3dp1 = res['dy3/dp1']
time = res['time']

To plot the trajectories using matplotlib, use the following commands:

plt.plot(time, dy1dp1, time, dy2dp1, time, dy3dp1)
plt.legend(('dy1/dp1', 'dy2/dp1', 'dy3/dp1'))
plt.show()

In Figure 7.10 the sensitivities are plotted.

This example will show how to use the JModelica.org's FMI-interface together with its simulation package, Assimulo. The FMU to be simulated is the full Robot from the Modelica standard library (3.1) where it is located in Mechanics.MultiBody.Examples.Systems.RobotR3. It consists of brakes, motors, gears and path planning. The model consists of 36 continuous states and around 700 algebraic variables together with 98 event functions and also a few thousand constants/parameters. The FMU was generated using Dymola 7.4.

This examples assumes that an FMU of the robot named:

Modelica_Mechanics_MultiBody_Examples_Systems_RobotR3_fullRobot.fmu

exists in the working folder.

Start by creating a Python script file and write or (copy paste) the command for importing the model object and the library used for plotting:

# Import the FMUModel class
from pyfmi import FMUModel

# Import the plotting library
import matplotlib.pyplot as plt

Next, we load the FMU into the model object:

robot = FMUModel('Modelica_Mechanics_MultiBody_Examples_Systems_RobotR3_fullRobot.fmu')

We are interested in simulating the Robot from time 0.0 to 1.8 using 1000 communication points and using tolerances specified in the FMU. This information is specified to the simulate method:

res = robot.simulate(start_time=0.0, final_time=1.8, options={'ncp':1000})

This preforms the simulation and the statistics will be printed in the prompt.

To retrieve data about a variable from the result data, access it as a dictionary with the name of the variable as key:

dq1  = res['der(mechanics.q[1])']
dq6  = res['der(mechanics.q[6])']
time = res['time']

Now we have loaded and retrieved the variables of interest. So lets plot them.

plt.plot(time,dq1,time,dq6)
plt.legend(['der(mechanics.q[1])','der(mechanics.q[6])'])
plt.xlabel('Time (s)')
plt.ylabel('Joint Velocity (rad/s)')
plt.title('Full Robot')
plt.show()

In Figure 7.12 the result is shown and in Figure 7.13 a comparison between Dymola and JModelica.org is plotted.