In the next sections, it will be shown how to use the JModelica.org platform for simulation of various JMUs and 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.
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 jmodelica.jmi import compile_jmu from jmodelica.jmi 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.
Simulation of ODEs in JModelica.org is currently limited to models written explicitly on the form:
There is also a limitation that events are currently not handled.
An example is the Van der Pol oscillator described in Section 2. In the code below it is shown again for convenience.
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;
We see here that the state derivatives are written on the left hand side while the equations are written on the right hand side.
Lets begin with the the example, copy and paste the Modelica code and save it into VDP.mo 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 jmodelica.jmi import compile_jmu from jmodelica.jmi import JMUModel # Import the plotting library import matplotlib.pyplot as plt
Next, we compile and load the model:
# Compile model
jmu_name = compile_jmu("VDP","VDP.mo")
# Load model
vdp = JMUModel(jmu_name)
Now we would like to change the default solver IDA to the ODE solver CVode and also change the relative tolerance to 0.0001. First we retrieve the default options and then change the option:
opts = vdp.simulate_options() #Retrieve the default options opts['solver'] = 'CVode' #Change the solver to CVode opts['CVode_options']['rtol'] = 0.0001 #Change the relative tolerance for CVode
Now lets provide the options to the simulate method and simulate the problem ten seconds:
res = vdp.simulate(final_time=10, options=opts)
The result is retrieved by accessing the res variable as a dictionary with the variable name as key:
x1 = res['x1'] x2 = res['x2'] t = res['time']
The variable trajectories are returned as numpy arrays and can be used for further analysis of the simulation result or for visualization:
plt.figure(1)
plt.plot(t, x1, t, x2)
plt.legend(('x1','x2'))
plt.title('Van der Pol oscillator.')
plt.ylabel('Angle (rad)')
plt.xlabel('Time (s)')
plt.show()
In Figure 7.5 the simulation result is shown.
The model which is to be simulated is an electric circuit. The model is depicted in Figure 7.6 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 jmodelica.jmi import compile_jmu from jmodelica.jmi 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.7.
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 jmodelica.jmi import compile_jmu, JMUModel
into the Python file. Compile and load the model:
# Define model file name and class name model_name = 'VDP' mofile = 'VDP.mo' # Compile model jmu_name = compile_jmu(model_name,mofile) # Load model vdp = JMUModel(jmu_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.8.
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 jmodelica.jmi import compile_jmu from jmodelica.jmi 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)
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.9 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 jmodelica.fmi 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.11 the result is shown and in Figure 7.12 a comparison between Dymola and JModelica.org is plotted.
