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()

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.)

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)

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 C.6, “General options for AssimuloAlg.” the general options for the `AssimuloAlg`

algorithm are described while in Table C.8, “Selection of solver arguments for IDA” a selection of the different solver arguments for the DAE solver IDA is shown. In Table C.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 C.6. General options for AssimuloAlg.**

Option | Default | Description |
---|---|---|

`solver` | `'IDA'` | Specifies the simulation method that is to be used. |

`ncp` | `0` | Number of communication points. If `ncp` is zero, the solver will return the internal steps taken. |

`initialize` | `True` | If 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_result` | `False` | Set 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_name` | Empty 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 C.7. Selection of solver arguments for CVode**

Option | Default | Description |
---|---|---|

`discr` | `'BDF'` | The discretization method. Can be either `'BDF'` or `'Adams'` |

`iter` | `'Newton'` | The iteration method. Can be either `'Newton` ' or `'FixedPoint'` . |

`maxord` | `5` | The maximum order used. Maximum for `'BDF'` is 5 while for the `'Adams'` method the maximum is 12 |

`maxh` | `Inf` | Maximum step-size. Positive float. |

`atol` | `1.0e-6` | Absolute 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. |

`rtol` | `1.0e-6` | Relative Tolerance. Positive float. |

**Table C.8. Selection of solver arguments for IDA**

Option | Default | Description |
---|---|---|

`maxord` | `5` | The maximum order used. Positive integer. |

`maxh` | `Inf` | Maximum step-size. Positive float. |

`atol` | `1.0e-6` | Absolute 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. |

`rtol` | `1.0e-6` | Relative Tolerance. Positive float. |

`suppress_alg` | `False` | Suppress the algebraic variables on the error test. Can be either `False` or `True` . |

`sensitivity` | `False` | If set to `True` , sensitivities for the states with respect to parameters set to free in the model will be calculated. |

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 C.9, “Result Object”.

**Table C.9. Result Object**

Option | Type | Description |
---|---|---|

`options` | Property | Gets the options object that was used during the simulation. |

`solver` | Property | Gets the solver that was used during the integration. |

`result_file` | Property | Gets the name of the generated result file. |

`is_variable(name)` | Method | Returns `True` if the given name is a time-varying variable. |

`data_matrix` | Property | Gets the raw data matrix. |

`is_negated(name)` | Method | Returns `True` if the given name is negated in the result matrix. |

`get_column(name)` | Method | Returns the column number in the data matrix which corresponds to the given variable. |

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

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 C.1, “A schematic picture of the quadruple tank process.” 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;

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 C.2, “Tank levels” the result of the tank levels are shown and in Figure C.3, “Input trajectories” the input signals are shown.

The model to be simulated in this example is an electric circuit. The model is depicted in Figure C.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`

.

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 C.5, “Simulation result”.

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 C.6, “Sensitivity results.” the sensitivities are plotted.