System Identification

System Identification (OMSysIdent) is part of the OpenModelica tool suite, but not bundled together with the main OpenModelica distribution and thus must be fetched separately from its project site.

OMSysIdent is a module for the parameter estimation for linear and nonlinear parametric dynamic models (wrapped as FMUs) on top of the OMSimulator API. It uses the Ceres solver (http://ceres-solver.org/) for the optimization task. The module provides a Python scripting API as well as an C API.

Note

Notice that this module was previously part of OMSimulator. It has been extracted out of the OMSimulator project and reorganized as a separate project in September 2020. As of 2020-10-07 the project is working on Linux but some more efforts are needed for migrating the Windows build and make the build and usage of the module more convenient.

Version: a65a0ed

Examples

There are examples in the testsuite which use the scripting API, as well as examples which directly use the C API.

Below is a basic example from the testsuite (HelloWorld_cs_Fit.py) which uses the Python scripting API. It determines the parameters for the following "hello world" style Modelica model:

model HelloWorld
  parameter Real a = -1;
  parameter Real x_start = 1;
  Real x(start=x_start, fixed=true);
equation
  der(x) = a*x;
end HelloWorld;

The goal is to estimate the value of the coefficent a and the initial value x_start of the state variable x. Instead of real measurements, the script simply uses simulation data generated from the HelloWorld examples as measurement data. The array data_time contains the time instants at which a sample is taken and the array data_x contains the value of x that corresponds to the respective time instant.

The estimation parameters are defined by calls to function simodel.addParameter(..) in which the name of the parameter and a first guess for the parameter's value is stated.

Listing 5 HelloWorld_cs_Fit.py

from OMSimulator import OMSimulator
from OMSysIdent import OMSysIdent
import numpy as np

oms = OMSimulator()

oms.setLogFile("HelloWorld_cs_Fit_py.log")
oms.setTempDirectory("./HelloWorld_cs_Fit_py/")
oms.newModel("HelloWorld_cs_Fit")
oms.addSystem("HelloWorld_cs_Fit.root", oms.system_wc)
# oms.setTolerance("HelloWorld_cs_Fit.root", 1e-5)

# add FMU
oms.addSubModel("HelloWorld_cs_Fit.root.HelloWorld", "../resources/HelloWorld.fmu")

# create simodel for model
simodel = OMSysIdent("HelloWorld_cs_Fit")
# simodel.describe()

# Data generated from simulating HelloWorld.mo for 1.0s with Euler and 0.1s step size
kNumSeries = 1
kNumObservations = 11
data_time = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1])
inputvars = []
measurementvars = ["root.HelloWorld.x"]
data_x = np.array([1, 0.9, 0.8100000000000001, 0.7290000000000001, 0.6561, 0.5904900000000001, 0.5314410000000001, 0.4782969000000001, 0.43046721, 0.387420489, 0.3486784401])

simodel.initialize(kNumSeries, data_time, inputvars, measurementvars)
# simodel.describe()

simodel.addParameter("root.HelloWorld.x_start", 0.5)
simodel.addParameter("root.HelloWorld.a", -0.5)
simodel.addMeasurement(0, "root.HelloWorld.x", data_x)
# simodel.describe()

simodel.setOptions_max_num_iterations(25)
simodel.solve("BriefReport")

status, state = simodel.getState()
# print('status: {0}; state: {1}').format(OMSysIdent.status_str(status), OMSysIdent.omsi_simodelstate_str(state))

status, startvalue1, estimatedvalue1 = simodel.getParameter("root.HelloWorld.a")
status, startvalue2, estimatedvalue2 = simodel.getParameter("root.HelloWorld.x_start")
# print('HelloWorld.a startvalue1: {0}; estimatedvalue1: {1}'.format(startvalue1, estimatedvalue1))
# print('HelloWorld.x_start startvalue2: {0}; estimatedvalue2: {1}'.format(startvalue2, estimatedvalue2))
is_OK1 = estimatedvalue1 > -1.1 and estimatedvalue1 < -0.9
is_OK2 = estimatedvalue2 > 0.9 and estimatedvalue2 < 1.1
print('HelloWorld.a estimation is OK: {0}'.format(is_OK1))
print('HelloWorld.x_start estimation is OK: {0}'.format(is_OK2))

# del simodel
oms.terminate("HelloWorld_cs_Fit")
oms.delete("HelloWorld_cs_Fit")

Running the script generates the following console output:

iter      cost      cost_change  |gradient|   |step|    tr_ratio  tr_radius  ls_iter  iter_time  total_time
 0  4.069192e-01    0.00e+00    2.20e+00   0.00e+00   0.00e+00  1.00e+04        0    7.91e-03    7.93e-03
 1  4.463938e-02    3.62e-01    4.35e-01   9.43e-01   8.91e-01  1.92e+04        1    7.36e-03    1.53e-02
 2  7.231043e-04    4.39e-02    5.16e-02   3.52e-01   9.85e-01  5.75e+04        1    7.26e-03    2.26e-02
 3  1.046555e-07    7.23e-04    4.74e-04   4.40e-02   1.00e+00  1.73e+05        1    7.31e-03    3.00e-02
 4  2.192358e-15    1.05e-07    5.77e-08   6.05e-04   1.00e+00  5.18e+05        1    7.15e-03    3.71e-02
 5  7.377320e-26    2.19e-15    2.05e-13   9.59e-08   1.00e+00  1.55e+06        1    7.42e-03    4.46e-02
Ceres Solver Report: Iterations: 6, Initial cost: 4.069192e-01, Final cost: 7.377320e-26, Termination: CONVERGENCE

=====================================
Total duration for parameter estimation: 44msec.
Result of parameter estimation (check 'Termination' status above whether solver converged):

HelloWorld_cs_Fit.root.HelloWorld.a(start=-0.5, *estimate*=-1)
HelloWorld_cs_Fit.root.HelloWorld.x_start(start=0.5, *estimate*=1)

=====================================
HelloWorld.a estimation is OK: True
HelloWorld.x_start estimation is OK: True
info:    Logging information has been saved to "HelloWorld_cs_Fit_py.log"

Python and C API

addInput

Add input values for external model inputs.

If there are several measurement series, all series need to be conducted with the same external inputs!

Python

Args:

var:

(str) Name of variable..

values:

(np.array) Array of input values for respective time instants in simodel.initialize().

Returns:

status:

(int) The C-API status code (oms_status_enu_t).

status = simodel.addInput(var, values)

C

oms_status_enu_t omsi_addInput(void* simodel, const char* var, const double* values, size_t nValues);

addMeasurement

Add measurement values for a fitting variable.

Python

Args:

iSeries:

(int) Index of measurement series.

var:

(str) Name of variable..

values:

(np.array) Array of measured values for respective time instants in simodel.initialize().

Returns:

status:

(int) The C-API status code (oms_status_enu_t).

status = simodel.addMeasurement(iSeries, var, values)

C

oms_status_enu_t omsi_addMeasurement(void* simodel, size_t iSeries, const char* var, const double* values, size_t nValues);

addParameter

Add parameter that should be estimated.

PYTHON

Args:

var:

(str) Name of parameter.

startvalue:

(float) Start value of parameter.

Returns:

status:

(int) The C-API status code (oms_status_enu_t).

status = simodel.addParameter(var, startvalue)

C

oms_status_enu_t omsi_addParameter(void* simodel, size_t iSeries, const char* var, const double* values, size_t nValues);

describe

Print summary of SysIdent model.

PYTHON

status = simodel.describe()

C

oms_status_enu_t omsi_describe(void* simodel);

freeSysIdentModel

Unloads a model.

PYTHON

Not available in Python. Related external C function called by class destructor.

C

void omsi_freeSysIdentModel(void* simodel);

getParameter

Get parameter that should be estimated.

PYTHON

Args:

var:

(str) Name of parameter.

Returns:

status:

(int) The C-API status code (oms_status_enu_t).

startvalue:

(float) Start value of parameter.

estimatedvalue:

(float) Estimated value of parameter.

status, startvalue, estimatedvalue = simodel.getParameter(var)

C

oms_status_enu_t omsi_getParameter(void* simodel, const char* var, double* startvalue, double* estimatedvalue);

getState

Get state of SysIdent model object.

PYTHON

Returns:

status:

(int) The C-API status code (oms_status_enu_t).

state:

(int) State of SysIdent model (omsi_simodelstate_t).

status, state = simodel.getState()

C

oms_status_enu_t omsi_getState(void* simodel, omsi_simodelstate_t* state);

initialize

This function initializes a given composite model. After this call, the model is in simulation mode.

PYTHON

Args:

nSeries:

(int) Number of measurement series.

time:

(numpy.array) Array of measurement/input time instants.

inputvars:

(list of str) List of names of input variables (empty list if none).

measurementvars:

(list of str) List of names of observed measurement variables.

Returns:

status:

(int) The C-API status code (oms_status_enu_t).

status = simodel.initalize(nSeries, time, inputvars, measurementvars)

C

oms_status_enu_t omsi_initialize(void* simodel, size_t nSeries, const double* time, size_t nTime, char const* const* inputvars, size_t nInputvars, char const* const* measurementvars, size_t nMeasurementvars);

newSysIdentModel

Creates an empty model for parameter estimation.

PYTHON

The corresponding Python function is the class constructor.

Args:

ident:

(str) Name of the model instance.

Returns:

simodel:

SysIdent model instance.

simodel = OMSysIdent(ident)

C

void* omsi_newSysIdentModel(const char* ident);

oms_status_str

Mapping of enum C-API status code (oms_status_enu_t) to string.

The C enum is reproduced below for convenience.

typedef enum {
  oms_status_ok,
  oms_status_warning,
  oms_status_discard,
  oms_status_error,
  oms_status_fatal,
  oms_status_pending
} oms_status_enu_t;

PYTHON

Args:

status:

(int) The C-API status code.

Returns:

status_str:

(str) String representation of status code.

The range of values of status corresponds to the C enum (by implicit conversion). This is a static Python method (@staticmethod).

status_str = oms_status_str(status)

C

Not available.

omsi_simodelstate_str

Mapping of enum C-API state code (omsi_simodelstate_t) to string.

The C enum is reproduced below for convenience.

typedef enum {
  omsi_simodelstate_constructed,    //!< After omsi_newSysIdentModel
  omsi_simodelstate_initialized,    //!< After omsi_initialize
  omsi_simodelstate_convergence,    //!< After omsi_solve if Ceres minimizer returned with ceres::TerminationType::CONVERGENCE
  omsi_simodelstate_no_convergence, //!< After omsi_solve if Ceres minimizer returned with ceres::TerminationType::NO_CONVERGENCE
  omsi_simodelstate_failure         //!< After omsi_solve if Ceres minimizer returned with ceres::TerminationType::FAILURE
} omsi_simodelstate_t;

PYTHON

Args:

state:

(int) State of SysIdent model.

Returns:

simodelstate_str:

(str) String representation of state code.

The range of values of state corresponds to the C enum (by implicit conversion). This is a static Python method (@staticmethod).

simodelstate_str = omsi_simodelstate_str(state)

C

Not available.

setOptions_max_num_iterations

Set Ceres solver option Solver::Options::max_num_iterations.

PYTHON

Args:

max_num_iterations:

(int) Maximum number of iterations for which the solver should run (default: 25).

Returns:

status:

(int) The C-API status code (oms_status_enu_t).

status = simodel.setOptions_max_num_iterations(max_num_iterations)

C

oms_status_enu_t omsi_setOptions_max_num_iterations(void* simodel, size_t max_num_iterations);

solve

Solve parameter estimation problem.

PYTHON

Args:

reporttype:

(str) Print report and progress information after call to Ceres solver. Supported report types: "", "BriefReport", "FullReport", where "" denotes no output.

Returns:

status:

(int) The C-API status code (oms_status_enu_t).

status = simodel.solve(reporttype)

C

oms_status_enu_t omsi_solve(void* simodel, const char* reporttype);