3. Simulation of JMUs

Simulation of JMUs in JModelica.org is performed via the simulate method of the JMU model object. The model object is called JMUModel and is located in the JModelica.org Python package pyjmi. JMUModel supports compiled models from JModelica.org which have the extension .jmu.

# Import JMUModel from pyjmi and load the JMU
from pyjmi import JMUModel
my_model = JMUModel('myJMU.jmu')

The simulation method in JMUModel is by default connected to the Assimulo simulation package and thus able to use its solvers. Continuing the short example from above, the following code will simulate the loaded JMU using default values and options:

res = my_model.simulate()

3.1. The simulate function

There are several parameters that can be set in the JMUModel.simulate function.

start_time

The time when the solver should start the integration.

final_time

The time when the solver should finish the integration.

input

Input signal for the simulation. (Further explained in below.)

algorithm

The algorithm that will be used for the simulation. Currently only a connection to Assimulo is supported and connected through the algorithm AssimuloAlg.

options

The options to be used in the algorithm. (Further explained in below.)

3.1.1. Input

The input defines the input trajectories to the model and should be a 2-tuple consisting of the name(s) of the input variables and the second argument should be either a data matrix or a function. If the argument is a data matrix it should contain a time vector as the first column and the second column should correspond to the first name in the first argument and so forth. If instead the second argument is a function it should be defined to take the time as input and return the number of inputs in the order defined by the first argument.

For example, consider that we have a model with an input variable u1 and that the model should be driven by a sinus wave as input. Also we are interested in the interval 0 to 10.

import numpy as N
t = N.linspace(0.,10.,100)            # Create one hundred evenly spaced points
u = N.sin(t)                          # Create the input vector
u_traj = N.transpose(N.vstack((t,u))) # Create the data matrix and transpose 
                                      # it to the correct form

The above code has created the data matrix that we are interested in giving to the model as input, we just need to connect the data to a specific input variable, u1:

input_object = ('u1', u_traj)

Now we are ready to simulate using the input and simulate 10 seconds.

res = model.simulate(final_time=10, input=input_object)

If we on the other hand would have two input variables, u1 and u2 the script would instead look like:

import numpy as N
t = N.linspace(0.,10.,100)                     # Create one hundred evenly spaced points
u1 = N.sin(t)                                  # Create the first input vector
u2 = N.cos(t)                                  # Create the second input vector
u_traj = N.transpose(N.vstack((t,u1,u2)))      # Create the data matrix and 
                                               # transpose it to the correct form
input_object = (['u1','u2'], u_traj)
res = model.simulate(final_time=10, input=input_object)

Note that the variables are now a List of variables.

If we were to do the same example using input functions instead, the code would look like for the single input case:

input_object = ('u1', N.sin)

and for the double input case:

def input_function(t):
    return N.array([N.sin(t),N.cos(t)])

input_object = (['u1','u2'],input_function)

3.1.2. Options for JMUModel

The options attribute are where options to the specified algorithm are stored and are preferably used together with:

opts = JMUModel.simulate_options()

which returns the default options for the default algorithm. Information about the available options can be viewed by typing help on the opts variable:

>>> help(opts)
   Options for simulation of a JMU model using the Assimulo simulation package.
   The Assimulo package contain both explicit solvers (CVode) for ODEs and 
   implicit solvers (IDA) for DAEs. The ODE solvers require that the problem
   is written on the form, ydot = f(t,y).
   
   ...

In Table 6.6, “General options for AssimuloAlg.” the general options for the AssimuloAlg algorithm are described while in Table 6.8, “Selection of solver arguments for IDA” a selection of the different solver arguments for the DAE solver IDA is shown. In Table 6.7, “Selection of solver arguments for CVode” a selection of solver arguments for the ODE solver CVode is shown. More information regarding the solver options can be found here, http://www.jmodelica.org/assimulo.

Table 6.6. General options for AssimuloAlg.

OptionDefaultDescription
solver'IDA'Specifies the simulation method that is to be used.
ncp0Number of communication points. If ncp is zero, the solver will return the internal steps taken.
initializeTrueIf set to True, an algorithm for initializing the differential equation is invoked, otherwise the differential equation is assumed to have consistent initial conditions.
write_scaled_resultFalseSet this parameter to True to write the result to file without taking scaling into account. If the value of scaled is False, then the variable scaling factors of the model are used to reproduced the unscaled variable values.
result_file_nameEmpty string (default generated file name will be used)Specifies the name of the file where the simulation result is written. Setting this option to an empty string results in a default file name that is based on the name of the model class.

Lets look at an example, consider that you want to simulate a JMU model using the solver CVode together with changing the discretization method (discr) from BDF to Adams:

...
opts = model.simulate_options()          # Retrieve the default options
opts['solver'] = 'CVode'                 # Change the solver from IDA to CVode
opts['CVode_options']['discr'] = 'Adams' # Change from using BDF to Adams
model.simulate(options=opts)             # Pass in the options to simulate and simulate

It should also be noted from the above example the options regarding a specific solver, say the tolerances for CVode or IDA, should be stored in a double dictionary where the first is named after the solver concatenated with _options:

opts['CVode_options']['atol'] = 1.0e-6 # Options specific for CVode
opts['IDA_options']['atol'] = 1.0e-6   # Options specific for IDA

For the general options, as changing the solver, they are accessed as a single dictionary:

opts['solver'] = 'CVode' # Changing the solver
opts['ncp'] = 1000       # Changing the number of communication points.

Table 6.7. Selection of solver arguments for CVode

OptionDefaultDescription
discr'BDF'The discretization method. Can be either 'BDF' or 'Adams'
iter'Newton'The iteration method. Can be either 'Newton' or 'FixedPoint'.
maxord5The maximum order used. Maximum for 'BDF' is 5 while for the 'Adams' method the maximum is 12
maxhInfMaximum step-size. Positive float.
atol1.0e-6Absolute Tolerance. Can be an array of floats where each value corresponds to the absolute tolerance for the corresponding variable. Can also be a single positive float.
rtol1.0e-6Relative Tolerance. Positive float.

Table 6.8. Selection of solver arguments for IDA

OptionDefaultDescription
maxord5The maximum order used. Positive integer.
maxhInfMaximum step-size. Positive float.
atol1.0e-6Absolute Tolerance. Can be an array of floats where each value corresponds to the absolute tolerance for the corresponding variable. Can also be a single positive float.
rtol1.0e-6Relative Tolerance. Positive float.
suppress_algFalseSuppress the algebraic variables on the error test. Can be either False or True.
sensitivityFalseIf set to True, sensitivities for the states with respect to parameters set to free in the model will be calculated.

3.1.3. Return argument

The return argument from the simulate method is an object derived from a common result object ResultBase in algorithm_drivers.py with a few extra convenience methods for retrieving the result of a variable. The result object can be accessed in the same way as a dictionary type in Python with the name of the variable as key.

res = model.simulate()
y = res['y']           # Return the result for the variable/parameter/constant y
dery = res['der(y)']   # Return the result for the variable/parameter/constant der(y)

This can be done for all the variables, parameters and constants defined in the model and is the preferred way of retrieving the result. There are however some more options available in the result object, see Table 6.9, “Result Object”.

Table 6.9. Result Object

OptionTypeDescription
optionsPropertyGets the options object that was used during the simulation.
solverPropertyGets the solver that was used during the integration.
result_filePropertyGets the name of the generated result file.
is_variable(name)MethodReturns True if the given name is a time-varying variable.
data_matrixPropertyGets the raw data matrix.
is_negated(name)MethodReturns True if the given name is negated in the result matrix.
get_column(name)MethodReturns the column number in the data matrix which corresponds to the given variable.

3.2. Examples

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

3.2.1. Simulation with inputs

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 also are the inputs to the model. The model is depicted in Figure 6.1, “A schematic picture of the quadruple tank process.” and in the code below the corresponding Modelica code is listed.

Figure 6.1. A schematic picture of the quadruple tank process.

A schematic picture of the quadruple tank process.

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;

Let's 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 Python/pyjmi/examples/files catalogue in the JModelica.org install folder. 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 has been adapted, let's 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 6.2, “Tank levels” the result of the tank levels are shown and in Figure 6.3, “Input trajectories” the input signals are shown.

Figure 6.2. Tank levels

Tank levels

Figure 6.3. Input trajectories

Input trajectories

3.2.2. Simulation of a discontinuous system

The model to be simulated in this example is an electric circuit. The model is depicted in Figure 6.4, “Electric Circuit” 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.

Figure 6.4. Electric Circuit

Electric Circuit

This example 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 in 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 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 6.5, “Simulation result”.

Figure 6.5. Simulation result

Simulation result

3.2.3. Simulation with sensitivities

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.

Let's begin with the the example. Copy and paste the Optimica code and save it into Robertson.mop, then 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)

Note that sensitivity computations are currently only supported for JMUModels. Now that the model is loaded, we have to change an option to activate the sensitivity calculations and also set 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 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 6.6, “Sensitivity results.” the sensitivities are plotted.

Figure 6.6. Sensitivity results.

Sensitivity results.