Changeset c35b1eb6 in OpenModelica


Ignore:
Timestamp:
2020-10-22T00:05:13+02:00 (3 years ago)
Author:
Bernhard Thiele <bernhard.thiele@…>
Children:
2041381
Parents:
df59bea9
git-author:
Bernhard Thiele <bernhard.thiele@…> (10/12/20 01:03:57)
git-committer:
Bernhard Thiele <bernhard.thiele@…> (10/22/20 00:05:13)
Message:

Migrate OMSysIdent docs into User's Guide

File:
1 edited

Legend:

Unmodified
Added
Removed
  • doc/UsersGuide/source/systemidentification.rst

    r1a316ea rc35b1eb6  
    22=====================
    33
    4 `System Identification (OMSysIdent) <https://github.com/OpenModelica/OMSysident>`_ is its own separate project and not bundled together with the main OpenModelica distribution.
     4`System Identification (OMSysIdent) <https://github.com/OpenModelica/OMSysident>`_
     5is part of the OpenModelica tool suite, but not bundled together with the main
     6OpenModelica distribution and thus must be fetched separately from its project site.
     7
     8*OMSysIdent* is a module for the parameter estimation for linear and nonlinear
     9parametric dynamic models (wrapped as FMUs) on top of the :doc:`omsimulator` API.
     10It uses the Ceres solver (http://ceres-solver.org/) for the optimization task.
     11The module provides a *Python scripting API* as well as an *C API*.
     12
     13.. note::
     14  Notice that this module was previously part of OMSimulator. It has been extracted
     15  out of the OMSimulator project and reorganized as a separate project in September 2020.
     16  As of 2020-10-07 the project is working on Linux but some more efforts are needed
     17  for migrating the Windows build and make the build and usage of the module
     18  more convenient.
     19
     20Version: `a65a0ed <https://github.com/OpenModelica/OMSysIdent/tree/a65a0edc3bdeebb1341fb3af8d3f100a4c86507a>`_
     21
     22Examples
     23########
     24
     25There are examples in the testsuite which use the scripting API, as well as
     26examples which directly use the C API.
     27
     28Below is a basic example from the testsuite (`HelloWorld_cs_Fit.py`) which
     29uses the Python scripting API. It determines the parameters for the following
     30"hello world" style Modelica model:
     31
     32.. code-block:: Modelica
     33
     34  model HelloWorld
     35    parameter Real a = -1;
     36    parameter Real x_start = 1;
     37    Real x(start=x_start, fixed=true);
     38  equation
     39    der(x) = a*x;
     40  end HelloWorld;
     41
     42The goal is to estimate the value of the coefficent `a` and the initial value
     43`x_start` of the state variable `x`. Instead of real measurements, the script
     44simply uses simulation data generated from the `HelloWorld` examples as
     45measurement data. The array `data_time` contains the time instants at which a
     46sample is taken and the array `data_x` contains the value of `x` that
     47corresponds to the respective time instant.
     48
     49The estimation parameters are defined by calls to function
     50`simodel.addParameter(..)` in which the name of the parameter and a first guess
     51for the parameter's value is stated.
     52
     53.. code-block:: python
     54  :caption: HelloWorld_cs_Fit.py
     55  :name: HelloWorld_cs_Fit-python
     56
     57  from OMSimulator import OMSimulator
     58  from OMSysIdent import OMSysIdent
     59  import numpy as np
     60
     61  oms = OMSimulator()
     62
     63  oms.setLogFile("HelloWorld_cs_Fit_py.log")
     64  oms.setTempDirectory("./HelloWorld_cs_Fit_py/")
     65  oms.newModel("HelloWorld_cs_Fit")
     66  oms.addSystem("HelloWorld_cs_Fit.root", oms.system_wc)
     67  # oms.setTolerance("HelloWorld_cs_Fit.root", 1e-5)
     68
     69  # add FMU
     70  oms.addSubModel("HelloWorld_cs_Fit.root.HelloWorld", "../resources/HelloWorld.fmu")
     71
     72  # create simodel for model
     73  simodel = OMSysIdent("HelloWorld_cs_Fit")
     74  # simodel.describe()
     75
     76  # Data generated from simulating HelloWorld.mo for 1.0s with Euler and 0.1s step size
     77  kNumSeries = 1
     78  kNumObservations = 11
     79  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])
     80  inputvars = []
     81  measurementvars = ["root.HelloWorld.x"]
     82  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])
     83
     84  simodel.initialize(kNumSeries, data_time, inputvars, measurementvars)
     85  # simodel.describe()
     86
     87  simodel.addParameter("root.HelloWorld.x_start", 0.5)
     88  simodel.addParameter("root.HelloWorld.a", -0.5)
     89  simodel.addMeasurement(0, "root.HelloWorld.x", data_x)
     90  # simodel.describe()
     91
     92  simodel.setOptions_max_num_iterations(25)
     93  simodel.solve("BriefReport")
     94
     95  status, state = simodel.getState()
     96  # print('status: {0}; state: {1}').format(OMSysIdent.status_str(status), OMSysIdent.omsi_simodelstate_str(state))
     97
     98  status, startvalue1, estimatedvalue1 = simodel.getParameter("root.HelloWorld.a")
     99  status, startvalue2, estimatedvalue2 = simodel.getParameter("root.HelloWorld.x_start")
     100  # print('HelloWorld.a startvalue1: {0}; estimatedvalue1: {1}'.format(startvalue1, estimatedvalue1))
     101  # print('HelloWorld.x_start startvalue2: {0}; estimatedvalue2: {1}'.format(startvalue2, estimatedvalue2))
     102  is_OK1 = estimatedvalue1 > -1.1 and estimatedvalue1 < -0.9
     103  is_OK2 = estimatedvalue2 > 0.9 and estimatedvalue2 < 1.1
     104  print('HelloWorld.a estimation is OK: {0}'.format(is_OK1))
     105  print('HelloWorld.x_start estimation is OK: {0}'.format(is_OK2))
     106
     107  # del simodel
     108  oms.terminate("HelloWorld_cs_Fit")
     109  oms.delete("HelloWorld_cs_Fit")
     110
     111Running the script generates the following console output:
     112
     113.. code-block:: none
     114
     115  iter      cost      cost_change  |gradient|   |step|    tr_ratio  tr_radius  ls_iter  iter_time  total_time
     116   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
     117   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
     118   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
     119   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
     120   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
     121   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
     122  Ceres Solver Report: Iterations: 6, Initial cost: 4.069192e-01, Final cost: 7.377320e-26, Termination: CONVERGENCE
     123
     124  =====================================
     125  Total duration for parameter estimation: 44msec.
     126  Result of parameter estimation (check 'Termination' status above whether solver converged):
     127
     128  HelloWorld_cs_Fit.root.HelloWorld.a(start=-0.5, *estimate*=-1)
     129  HelloWorld_cs_Fit.root.HelloWorld.x_start(start=0.5, *estimate*=1)
     130
     131  =====================================
     132  HelloWorld.a estimation is OK: True
     133  HelloWorld.x_start estimation is OK: True
     134  info:    Logging information has been saved to "HelloWorld_cs_Fit_py.log"
     135
     136Python and C API
     137################
     138
     139addInput
     140--------
     141
     142Add input values for external model inputs.
     143
     144If there are several measurement series, all series need to be conducted
     145with the same external inputs!
     146
     147
     148Python
     149^^^^^^
     150
     151Args:
     152  :var: (str) Name of variable..
     153  :values: (np.array) Array of input values for respective time instants in `simodel.initialize()`.
     154
     155Returns:
     156  :status: (int) The C-API status code (`oms_status_enu_t`).
     157
     158.. code-block:: python
     159
     160  status = simodel.addInput(var, values)
     161
     162
     163C
     164^
     165
     166.. code-block:: c
     167
     168  oms_status_enu_t omsi_addInput(void* simodel, const char* var, const double* values, size_t nValues);
     169
     170
     171addMeasurement
     172--------------
     173
     174Add measurement values for a fitting variable.
     175
     176Python
     177^^^^^^
     178
     179Args:
     180  :iSeries: (int) Index of measurement series.
     181  :var: (str) Name of variable..
     182  :values: (np.array) Array of measured values for respective time instants in `simodel.initialize()`.
     183
     184Returns:
     185  :status: (int) The C-API status code (`oms_status_enu_t`).
     186
     187.. code-block:: python
     188
     189  status = simodel.addMeasurement(iSeries, var, values)
     190
     191C
     192^
     193
     194.. code-block:: c
     195
     196  oms_status_enu_t omsi_addMeasurement(void* simodel, size_t iSeries, const char* var, const double* values, size_t nValues);
     197
     198
     199addParameter
     200------------
     201
     202Add parameter that should be estimated.
     203
     204PYTHON
     205^^^^^^
     206
     207Args:
     208  :var: (str) Name of parameter.
     209  :startvalue: (float) Start value of parameter.
     210
     211Returns:
     212  :status: (int) The C-API status code (`oms_status_enu_t`).
     213
     214.. code-block:: python
     215
     216  status = simodel.addParameter(var, startvalue)
     217
     218C
     219^
     220
     221.. code-block:: c
     222
     223  oms_status_enu_t omsi_addParameter(void* simodel, size_t iSeries, const char* var, const double* values, size_t nValues);
     224
     225
     226describe
     227--------
     228
     229Print summary of SysIdent model.
     230
     231PYTHON
     232^^^^^^
     233
     234.. code-block:: python
     235
     236  status = simodel.describe()
     237
     238C
     239^
     240
     241.. code-block:: c
     242
     243  oms_status_enu_t omsi_describe(void* simodel);
     244
     245
     246freeSysIdentModel
     247-----------------
     248
     249Unloads a model.
     250
     251PYTHON
     252^^^^^^
     253
     254Not available in Python. Related external C function called by class destructor.
     255
     256
     257C
     258^
     259
     260.. code-block:: c
     261
     262  void omsi_freeSysIdentModel(void* simodel);
     263
     264
     265getParameter
     266------------
     267
     268Get parameter that should be estimated.
     269
     270PYTHON
     271^^^^^^
     272
     273Args:
     274  :var: (str) Name of parameter.
     275
     276Returns:
     277  :status: (int) The C-API status code (`oms_status_enu_t`).
     278  :startvalue: (float) Start value of parameter.
     279  :estimatedvalue: (float) Estimated value of parameter.
     280
     281.. code-block:: python
     282
     283  status, startvalue, estimatedvalue = simodel.getParameter(var)
     284
     285C
     286^
     287
     288.. code-block:: c
     289
     290  oms_status_enu_t omsi_getParameter(void* simodel, const char* var, double* startvalue, double* estimatedvalue);
     291
     292
     293getState
     294--------
     295
     296Get state of SysIdent model object.
     297
     298PYTHON
     299^^^^^^
     300
     301Returns:
     302  :status: (int) The C-API status code (`oms_status_enu_t`).
     303  :state: (int) State of SysIdent model (`omsi_simodelstate_t`).
     304
     305.. code-block:: python
     306
     307  status, state = simodel.getState()
     308
     309C
     310^
     311
     312.. code-block:: c
     313
     314  oms_status_enu_t omsi_getState(void* simodel, omsi_simodelstate_t* state);
     315
     316
     317initialize
     318----------
     319
     320This function initializes a given composite model. After this call, the model is in simulation mode.
     321
     322PYTHON
     323^^^^^^
     324
     325Args:
     326  :nSeries: (int) Number of measurement series.
     327  :time: (numpy.array) Array of measurement/input time instants.
     328  :inputvars: (list of str) List of names of input variables (empty list if none).
     329  :measurementvars: (list of str) List of names of observed measurement variables.
     330
     331Returns:
     332  :status: (int) The C-API status code (`oms_status_enu_t`).
     333
     334.. code-block:: python
     335
     336  status = simodel.initalize(nSeries, time, inputvars, measurementvars)
     337
     338C
     339^
     340
     341.. code-block:: c
     342
     343  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);
     344
     345
     346newSysIdentModel
     347----------------
     348
     349Creates an empty model for parameter estimation.
     350
     351PYTHON
     352^^^^^^
     353
     354The corresponding Python function is the class constructor.
     355
     356Args:
     357  :ident: (str) Name of the model instance.
     358
     359Returns:
     360  :simodel: SysIdent model instance.
     361
     362.. code-block:: python
     363
     364  simodel = OMSysIdent(ident)
     365
     366C
     367^
     368
     369.. code-block:: c
     370
     371  void* omsi_newSysIdentModel(const char* ident);
     372
     373
     374setOptions_max_num_iterations
     375-----------------------------
     376
     377Set Ceres solver option `Solver::Options::max_num_iterations`.
     378
     379PYTHON
     380^^^^^^
     381
     382Args:
     383  :max_num_iterations: (int) Maximum number of iterations for which the solver should run (default: 25).
     384
     385Returns:
     386  :status: (int) The C-API status code (`oms_status_enu_t`).
     387
     388.. code-block:: python
     389
     390  status = simodel.setOptions_max_num_iterations(max_num_iterations)
     391
     392C
     393^
     394
     395.. code-block:: c
     396
     397  oms_status_enu_t omsi_setOptions_max_num_iterations(void* simodel, size_t max_num_iterations);
     398
     399
     400solve
     401-----
     402
     403Solve parameter estimation problem.
     404
     405PYTHON
     406^^^^^^
     407
     408Args:
     409  :reporttype: (str) Print report and progress information after call to Ceres solver.
     410               Supported report types: `"", "BriefReport", "FullReport"`, where `""` denotes no output.
     411
     412Returns:
     413  :status: (int) The C-API status code (`oms_status_enu_t`).
     414
     415.. code-block:: python
     416
     417  status = simodel.solve(reporttype)
     418
     419C
     420^
     421
     422.. code-block:: c
     423
     424  oms_status_enu_t omsi_solve(void* simodel, const char* reporttype);
Note: See TracChangeset for help on using the changeset viewer.