Sensitivity Analysis using PyFMI - FMU in for-loop

Question

Main goal

Sensitivity analysis of a district heating network.

Approach

1. Modelica model of the system (in Dymola) using the AixLib and BuildingSystem libraries

2. Export model as FMU co-simulation

3. Use SALib (sensitivity analysis python library) to define the samples (parameter sweep)

4. Use PyFMI to run the model in a for-loop in Python for all the individual samples (and parallelize the for loop maybe using JobLib to perfome the simulation on multiple processors)

5. SALib to perform a variance-based sensitivity analyses (http://salib.readthedocs.io/en/latest/basics.html#an-example)

First step

Simple modelica model of the Ishigami function (not time dependent). This function is often used to test sensitivity analysis methods (https://www.sfu.ca/~ssurjano/ishigami.html).

The python code (including loading the FMU with PyFMI and the parameter sweep) works fine.

The problem

After a certain amount of simulation we get an error. The error output looks not always the same. Sometimes we get

FMUException: Error loading the binary. Could not load the DLL: Eine DLL-Initialisierungsroutine ist fehlgeschlagen.

Translation: A DLL-Initilisation routine is failed.

And sometimes we get:

FMUException: Error loading the binary. Could not load the DLL: Für diesen Befehl ist nicht genügend Speicher verfügbar.

Translation: There is not enough memory available for this command.

The error occurs after around 650 simulation runs. This is not dependent of if the simulations are performed in smaller loop-blocks which are re-run one after another or if one single for loop runs through all the simulations. By restarting the python console/process, new simulations can be run again.

Working environment:

Windows 10, Python 2.7, PyFMI installed using pip (not JModelica), Python coding on Jupyther notebook (on Mozilla Firefox)

We have only basic knowledge of python and PyFMI and are really struggling with this error!

Attachment

Below you can find

• Modelica model used to export co-simulation FMU from Dymola (using CVode)

• Python code as py file

• Output scatter plot of the python code.

I also made a post on the JModelica Forum, where you can download the files directly (FMU, Jupyter notebook, etc.): http://www.jmodelica.org/27925

Modelica model

model IshigamiFunction

  final parameter Real a = 7;

  final parameter Real b = 0.05;

  parameter Real x1 = 1;

  parameter Real x2 = 1;

  parameter Real x3 = 1;

  Real f;

equation

  f = sin(x1) + a * sin(x2)^2 + b * x3^4 * sin(x1);

end IshigamiFunction;

Python code

import numpy as np
import pylab as pl
from pyfmi import load_fmu
from SALib.sample import saltelli
from SALib.analyze import sobol
from ipywidgets import FloatProgress
from IPython.display import display


n = 100

problem = {
    'num_vars': 3,
    'names': ['x1', 'x2', 'x3'],
    'bounds': [[-np.pi, np.pi],
               [-np.pi, np.pi],
               [-np.pi, np.pi]]
}

param_values = saltelli.sample(problem, n)

fmu = 'Model\IshigamiFunction\IshigamiFunction.fmu'
n_sim = param_values.shape[0]

# Progress bar
f = FloatProgress(min = 0, max = n_sim, description='Progress:')
display(f)

# Numpy array to save results
y = np.zeros([param_values.shape[0]])
x1 = np.zeros([param_values.shape[0]])
x2 = np.zeros([param_values.shape[0]])
x3 = np.zeros([param_values.shape[0]])

for i, X in enumerate(param_values):
    model = load_fmu(fmu)  
    model.set(problem['names'], X)
    res = model.simulate(final_time = 1)
    y[i] = res['f'][-1]
    x1[i] = res['x1'][-1]
    x2[i] = res['x2'][-1]
    x3[i] = res['x3'][-1]
    f.value += 1


# Scatter plots
fig = pl.figure(figsize=(20, 5))
pl.clf()

pl.subplot(1,3,1)
pl.plot(x1, y, 'or')
pl.ylabel('x1')
pl.xlabel('f')

pl.subplot(1,3,2)
pl.plot(x2, y, 'ob')
pl.ylabel('x2')
pl.xlabel('f')

pl.subplot(1,3,3)
pl.plot(x3, y, 'og')
pl.ylabel('x3')
pl.xlabel('f')

pl.suptitle('Scatter plots')
pl.show()

# Sensitivity analysis
Si = sobol.analyze(problem, y, print_to_console=True)

Output plot from python script

Update

I made some more tests, and this is what I found:

Depending on if the FMU is exported from Dymola or from JModelica the behavior is different:

Using an FMU exported from Dymola:

• Taking the load_fmu line out of the for-loop seems to work
• Even with the load_fmu not in the for-loop there are sometimes crashes
• Adding a new line model.reset() before the model.set(...) command seems to work fine
• The results are different when simulated with or without model.reset() -> Why??
• model.instantiate() instead of model.reset() -> doesn't work. The memory usage in the task manager goes up to around 350 MB and then

  FMUException: Failed to instantiate the model. See the log for possibly more information.

The log file with log_level=4:

FMIL: module = FMILIB, log level = 4: XML specifies FMI standard version 2.0
FMIL: module = FMILIB, log level = 4: Loading 'win32' binary with 'default' platform types
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiInstantiateModel completed
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiInstantiateSlave
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiInstantiateModel completed
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiInstantiateSlave
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetReal: x1 = -1.76101
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetReal: x2 = -2.53414
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetReal: x3 = 0.116583
FMIL: module = Model, log level = 4: [][FMU status:OK] fmi2SetupExperiment: startTime is set to 0
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiEnterSlaveInitializationMode...
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiEnterSlaveInitializationMode completed
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiExitSlaveInitializationMode...
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiExitSlaveInitializationMode completed
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetReal: x1 = -1.76101
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetReal: x2 = -2.53414
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetReal: x3 = 0.116583
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetReal: a = 7
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetReal: b = 0.05
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetReal: f = 1.29856
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetDerivatives
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetDerivatives
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetDerivatives
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetDerivatives
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetDerivatives
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetDerivatives
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetDerivatives
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetDerivatives
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.002
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.004
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.006
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiGetDerivatives
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.008
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.01
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.012
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.014
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.016
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.018
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.02
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep

...

FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.99
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.992
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.994
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.996
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 0.998
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiSetTime to 1
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiDoStep
FMIL: module = Model, log level = 1: [][FMU status:Fatal] The license file was not found. Use the environment variable "DYMOLA_RUNTIME_LICENSE" to specify your Dymola license file.

FMIL: module = Model, log level = 1: [][FMU status:Fatal] Instantiation failed
FMIL: module = Model, log level = 4: [][FMU status:OK] fmiFreeModelInstance

Using an FMU exported from JModelica:

• Works fine even if the load_fmu is within the for-loop (but slower)
• This experience does not correspond the example given within the JModelica documentation in chapter 5.4.2 (http://www.jmodelica.org/api-docs/usersguide/2.1/ch05s04.html#d0e1854) where the load_fmu command is given within the for-loop
• The command model.reset() or model.instatiate() is required within the for-loop (contrary to Dymola FMU) -> Why??

My questions:

What is the right why to perform a loop, which simulates a FMU model many times with differen parameters?

What is the difference between using model.reset(), model.instatiate() or none of them?

Attachment

Here's a plot showing the difference between a for-loop with model.reset() and without it.

An FMU exported from JModelica (doesn't need any license) can be downloaded here: http://www.jmodelica.org/27925#comment-6668

Answer 1

The right way for a Dymola FMU (and probably the same for FMUs from other vendors) would be call fmi/fmi2Instantiate outside the for loop. These functions will allocate memory and perform license check if the FMU is exported without binary export license. By calling fmiResetSlave/fmi2Reset you can reset the FMU to the instantiated state without new memory allocation.

• fmiInstantiateSlave/fmi2Instantiate

  creates an FMU instance that can be used for simulation, multiple calls will create multiple instances, each which would require new memory allocation and proper deletion.

• fmiReset

  resets your instance to the state after instantiate and before fmiInitializeSlave/fmi2Intialize is called. This is faster requires no new dynamic memory allocation and should be used in your case.

In addition the license check of Dymola FMU's exported without binary export may leak memory in older Dymola versions when calling fmiFreeSalveInstance/fmi2FreeInstance. This is in most cases not a problem as you normally terminate your program when you terminate your FMU. By instantiating your FMU inside the for loop this becomes serious and your memory will finally end. Fix pack should be available if you contact Dymola support.

Answer 2

It looks like a memory issue to me. Can you observe the allocated memory during your run in the Win Task Manager? Btw, your FMU (from your cross post) requires a DYMOLA_RUNTIME_LICENSE, which limits the reproduction to Dymola users only.