Dymola Simulation works, FMPY not

[weberscode]

Good day to all,

I have the following problem: the simulation in Dymola with CVODE works, with DASSL it does not. The model is very stiff. I would like to do a parameter variation in python and other things based on it. Unfortunately I can not simulate the CoSimulationFMU generated by Dymola with CVODE solvers. What am I doing wrong? Do I need to get the model running in Dymola with Dassl first? I suspect so, since the fmu simulation in python and the model_01_Implementations_BS.zip Dassl simulation in Dymola fail at the same point. I am very grateful for any answer that will advance our/my research. KR Simon

[t-sommer]

Have you tried to export the FMU as "Co-Simulation using Cvode"?

[weberscode]

Thank you for the quick answer Mr. Sommer,

I have exported the FMU only as "Co-Simulation using CVODe".

My problem is that I can simulate the simulation in Dymola using CVODE, then export it and then only use it with python - fmpy in a very limited way.

Some simulations do show fmi2getDerivatives: dsblock_failed, QiErr = 1, but then runs through.

Some simulations show something and abort.

So it is not possible for me to carry out the parameter variation. What could I change? I have attached the code.

# =============================================================================
# LIBS
# =============================================================================
import numpy as np
import pandas as pd
import os
import sys
import fmpy
from fmpy.fmi2 import FMU2Slave
from fmpy import read_model_description
from matplotlib import pyplot as plt
from scipy.optimize import minimize
from scipy.optimize import minimize_scalar
from scipy.optimize import Bounds

# =============================================================================
# Functions
# =============================================================================
 

    
# =============================================================================
# INIT parameters
# =============================================================================

filename = 'model_01_Implementations_BS_sch.fmu' # Filename

# check file availability
fmu_path = 'FMU/'
dymola_work_path = 'C:/Users/siosw/Documents/01_Promotion/02_C06/Dymola2/work/'
fmu_filename = fmu_path + filename
dymola_work_path_filename = dymola_work_path + filename
if (os.path.exists(fmu_filename) and os.path.exists(dymola_work_path_filename) ) or os.path.exists(dymola_work_path_filename):
    os.replace(dymola_work_path_filename,fmu_filename)
elif os.path.exists(fmu_filename):
    pass
else:
    print("error: No fmu file found with specified name", file=sys.stderr)

# read the model description
model_description = read_model_description(fmu_filename)

# collect the value references for the variables to read / write
vrs = {}
for variable in model_description.modelVariables:
    vrs[variable.name] = variable.valueReference

# extract the FMU
unzipdir = fmpy.extract(fmu_filename)

fmu_args = {'guid': model_description.guid,
            'modelIdentifier': model_description.coSimulation.modelIdentifier,
            'unzipDirectory': unzipdir}

# get the value references for the start and output values
vrs_swi_real = [vrs['partialC.det_u_ADCH_c.sel_phase1.sel_phase_sd_sch_2_1.t23'],vrs['partialC.det_u_ADCH_c.sel_phase1.sel_phase_sd_sch_2_1.h2'],vrs['partialC.det_u_ADCH_c.sel_phase1.sel_phase_sd_sch_2_1.h3']]
vrs_sta_real_init = [vrs['init.T_z_0'],vrs['init.T_e_0'],vrs['init.T_c_0'],vrs['init.T_a_0'],vrs['init.X_0']]
vrs_sta_real_end = [vrs['partialP.A.bbus.T_Z'],vrs['partialP.A.hbus.T_ADCH[2]'],vrs['partialP.A.hbus.T_ADCH[3]'],vrs['partialP.A.hbus.T_ADCH[1]'],vrs['partialP.A.hbus.X']]
vrs_par_int = [vrs['site.or_a[1]'],vrs['site.or_c[1]']]
vrs_par_real = [vrs['site.f_a[1]'],vrs['site.f_c[1]'],vrs['site.f_e'],vrs['site.azi']]
vrs_res = [vrs['partialP.A.abus.E_tot']]

fmu = FMU2Slave(**fmu_args)
#dump(fmu)  # get information

fmu.instantiate(visible=True)

"""
FUNCTIONS
"""
def get_periodicity(args):
    vrs_sta_real_init, val_sta_real_old, vrs_swi_real, val_swi_real = args
    n = 0
    nmax = 100
    tol_T = 1e-2
    tol_X = 1e-4
    tol = np.array([tol_T,tol_T,tol_T,tol_T,tol_X])
    residual = np.ones(5)/tol
    print("----------------------")
    print("get periodic beginnt")
    print(residual.round(decimals=2))  
    print_val = val_sta_real_old-np.array([T0,T0,T0,T0,0])
    print_val[0:4] = print_val[0:4].round(decimals=2)
    print_val[4] = print_val[4].round(decimals=5)
    print(print_val)  
    
    while any(residual>np.ones(5)) and n<nmax:
        fmu.reset()
 
        fmu.setupExperiment(tolerance=1e-8)
    
        # set the start values
        fmu.setReal(vr=vrs_sta_real_init, value=val_sta_real_old)
        # set the switching time t23
        fmu.setReal(vr=vrs_swi_real, value=val_swi_real)
        
        fmu.enterInitializationMode()
        fmu.exitInitializationMode()
    
        start = 0.0
        end = 86400
    
        # simulation loop
        #print('Im here')
        fmu.doStep(currentCommunicationPoint=start, communicationStepSize=end)
    
        # get the start values for the next loop
        val_sta_real_new = np.array(fmu.getReal(vrs_sta_real_end))
        val_sta_real_new[1:4] = val_sta_real_new[1:4]+T0
        residual = abs(val_sta_real_new-val_sta_real_old)/tol
        #val_sta_real_old = (val_sta_real_old + val_sta_real_new)/2
        val_sta_real_old = val_sta_real_new
        n = n+1
        print('while loop nr' + str(n))
    
    print(residual.round(decimals=2))  
    print_val = val_sta_real_old-np.array([T0,T0,T0,T0,0])
    print_val[0:4] = print_val[0:4].round(decimals=2)
    print_val[4] = print_val[4].round(decimals=5)
    print(print_val)  
    print("get periodic abgeschlossen")
    print("----------------------")
    return val_sta_real_old

def get_min_energy_at_swi_1(x,*args):
    vrs_sta_real_init, val_sta_real_old, vrs_swi_real, val_swi_real = args
    #vrs_sta_real_init, val_sta_real_old, vrs_swi_real = args
    
    fmu.reset()

    fmu.setupExperiment(tolerance=1e-8)

    # set the initial state values
    fmu.setReal(vr=vrs_sta_real_init, value=val_sta_real_old)
    # set the switching time t23
    val_swi_real[0] = x
    fmu.setReal(vr=vrs_swi_real, value=val_swi_real)

    fmu.enterInitializationMode()
    fmu.exitInitializationMode()

    start = 0.0
    end = 86400

    # simulation loop
    #print('Im here')
    fmu.doStep(currentCommunicationPoint=start, communicationStepSize=end)    
    res = np.array(fmu.getReal(vrs_res))
    return res*1000

def get_min_energy_at_swi_2(x,*args):
    vrs_sta_real_init, val_sta_real_old, vrs_swi_real, val_swi_real = args
    #vrs_sta_real_init, val_sta_real_old, vrs_swi_real = args
    
    fmu.reset()

    fmu.setupExperiment(tolerance=1e-8)

    # set the initial state values
    fmu.setReal(vr=vrs_sta_real_init, value=val_sta_real_old)
    # set the switching time t23
    val_swi_real[1] = x
    fmu.setReal(vr=vrs_swi_real, value=val_swi_real)

    fmu.enterInitializationMode()
    fmu.exitInitializationMode()

    start = 0.0
    end = 86400

    # simulation loop
    #print('Im here')
    fmu.doStep(currentCommunicationPoint=start, communicationStepSize=end)    
    res = np.array(fmu.getReal(vrs_res))
    return res*1000

def get_min_energy_at_swi_3(x,*args):
    
    
    vrs_sta_real_init, val_sta_real_old, vrs_swi_real, val_swi_real = args
    #vrs_sta_real_init, val_sta_real_old, vrs_swi_real = args
    
    fmu.reset()

    fmu.setupExperiment(tolerance=1e-8)

    # set the initial state values
    fmu.setReal(vr=vrs_sta_real_init, value=val_sta_real_old)
    # set the switching time t23
    val_swi_real[2] = x
    fmu.setReal(vr=vrs_swi_real, value=val_swi_real)

    fmu.enterInitializationMode()
    fmu.exitInitializationMode()

    start = 0.0
    end = 86400

    # simulation loop
    #print('Im here')
    fmu.doStep(currentCommunicationPoint=start, communicationStepSize=end)    
    res = np.array(fmu.getReal(vrs_res))
    return res*1000

def get_min_steady_state_energy(args):
    for i in range(3):  
        print("for loop " + str(i))
        
        val_sta_real_old = get_periodicity(args)
        args=(vrs_sta_real_init, val_sta_real_old, vrs_swi_real,val_swi_real) # acutalize initial states
        
        
        # optimize Etot under t23
        print(val_swi_real)
        
        res = minimize(get_min_energy_at_swi_1, val_swi_real[0], method='trust-constr',args=args,tol=0.01,bounds=bnds_sw1)
        val_swi_real[0]=res.x
        args=(vrs_sta_real_init, val_sta_real_old, vrs_swi_real,val_swi_real)
        print(res)
        print(val_swi_real)
        
        res = minimize(get_min_energy_at_swi_2, val_swi_real[1], method='trust-constr',args=args,tol=0.01,bounds=bnds_sw2)
        val_swi_real[1]=res.x
        args=(vrs_sta_real_init, val_sta_real_old, vrs_swi_real,val_swi_real)
        print(res)
        print(val_swi_real)
        
        res = minimize(get_min_energy_at_swi_3, val_swi_real[2], method='trust-constr',args=args,tol=0.01,bounds=bnds_sw3)
        val_swi_real[2]=res.x
        args=(vrs_sta_real_init, val_sta_real_old, vrs_swi_real,val_swi_real)
        print(res)
        print(val_swi_real)
    
    return args

    
# iterate until stationary criterias are fullfilled
#sta_names = ['T_z','T_e','T_c','T_a','X']
#par_names = ['or_a','or_c','f_a','f_c','f_e','azi','E_tot']
#par_val_init = np.array([1,4,0.5,0.5,0.5,-20]) # initial parameter point
T0 = 273.15
val_sta_real_old = np.array([24.02+T0,24.02+T0,25.7+T0,49+T0,0.2641]) # initial state point
val_swi_real = np.array([13.117,3.526,2.144]) # initial switch times
bnds_sw1 = Bounds(9.0, 14.0)
bnds_sw2 = Bounds(1.0, 5.0)
bnds_sw3 = Bounds(1.0, 5.0)
         
args=(vrs_sta_real_init, val_sta_real_old, vrs_swi_real,val_swi_real)


args = get_min_steady_state_energy(args)


vrs_sta_real_init, val_sta_real_old, vrs_swi_real, val_swi_real = args


fmu.terminate()

# call the FMI API directly to avoid unloading the share library
fmu.fmi2FreeInstance(fmu.component)

    
# sta_names = ['T_z','T_e','T_c','T_a','X']
# par_names = ['or_a','or_c','f_a','f_c','f_e','azi','E_tot']
# par_val_init = np.array([1,4,0.5,0.5,0.5,-20]) # initial parameter point
# sta_val_init = np.array([23,23,20,20,0.335]) # initial state point
# res_val_init = np.zeros(1)
# df_stor = pd.DataFrame(columns=par_names)
# df_stor.loc[0,:] = np.concatenate((par_val_init,res_val_init))



# =============================================================================
# get steady state initial state values to parameter point
# ============================================================================= 

[t-sommer]

@karlwernersson, do you have an idea what might cause this behavior?

[weberscode]

Currently I Think it is on the model side. I got around the problem by try: except:. In except case: i adjust the initial point, restart and it works quite well for me.

[KarlWernersson]

Sorry for the late reply The model seems very sensitive as it fails with dassl but succeds with cvode. So my genral sudgestion would be to inveistage the model and why it fails in Dymola to see if it can be made more robust.

Then the reason for failure when simulating the cs fmu in fmpy vs dymola.

if no input. tolerance might be set differently.

if input when input is given and size of doStep can give variation that could affect a very sensitive model.