Changeset c35b1eb6 in OpenModelica
- Timestamp:
- 2020-10-22T00:05:13+02:00 (3 years ago)
- 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)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
doc/UsersGuide/source/systemidentification.rst
r1a316ea rc35b1eb6 2 2 ===================== 3 3 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>`_ 5 is part of the OpenModelica tool suite, but not bundled together with the main 6 OpenModelica 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 9 parametric dynamic models (wrapped as FMUs) on top of the :doc:`omsimulator` API. 10 It uses the Ceres solver (http://ceres-solver.org/) for the optimization task. 11 The 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 20 Version: `a65a0ed <https://github.com/OpenModelica/OMSysIdent/tree/a65a0edc3bdeebb1341fb3af8d3f100a4c86507a>`_ 21 22 Examples 23 ######## 24 25 There are examples in the testsuite which use the scripting API, as well as 26 examples which directly use the C API. 27 28 Below is a basic example from the testsuite (`HelloWorld_cs_Fit.py`) which 29 uses 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 42 The 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 44 simply uses simulation data generated from the `HelloWorld` examples as 45 measurement data. The array `data_time` contains the time instants at which a 46 sample is taken and the array `data_x` contains the value of `x` that 47 corresponds to the respective time instant. 48 49 The estimation parameters are defined by calls to function 50 `simodel.addParameter(..)` in which the name of the parameter and a first guess 51 for 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 111 Running 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 136 Python and C API 137 ################ 138 139 addInput 140 -------- 141 142 Add input values for external model inputs. 143 144 If there are several measurement series, all series need to be conducted 145 with the same external inputs! 146 147 148 Python 149 ^^^^^^ 150 151 Args: 152 :var: (str) Name of variable.. 153 :values: (np.array) Array of input values for respective time instants in `simodel.initialize()`. 154 155 Returns: 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 163 C 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 171 addMeasurement 172 -------------- 173 174 Add measurement values for a fitting variable. 175 176 Python 177 ^^^^^^ 178 179 Args: 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 184 Returns: 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 191 C 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 199 addParameter 200 ------------ 201 202 Add parameter that should be estimated. 203 204 PYTHON 205 ^^^^^^ 206 207 Args: 208 :var: (str) Name of parameter. 209 :startvalue: (float) Start value of parameter. 210 211 Returns: 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 218 C 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 226 describe 227 -------- 228 229 Print summary of SysIdent model. 230 231 PYTHON 232 ^^^^^^ 233 234 .. code-block:: python 235 236 status = simodel.describe() 237 238 C 239 ^ 240 241 .. code-block:: c 242 243 oms_status_enu_t omsi_describe(void* simodel); 244 245 246 freeSysIdentModel 247 ----------------- 248 249 Unloads a model. 250 251 PYTHON 252 ^^^^^^ 253 254 Not available in Python. Related external C function called by class destructor. 255 256 257 C 258 ^ 259 260 .. code-block:: c 261 262 void omsi_freeSysIdentModel(void* simodel); 263 264 265 getParameter 266 ------------ 267 268 Get parameter that should be estimated. 269 270 PYTHON 271 ^^^^^^ 272 273 Args: 274 :var: (str) Name of parameter. 275 276 Returns: 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 285 C 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 293 getState 294 -------- 295 296 Get state of SysIdent model object. 297 298 PYTHON 299 ^^^^^^ 300 301 Returns: 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 309 C 310 ^ 311 312 .. code-block:: c 313 314 oms_status_enu_t omsi_getState(void* simodel, omsi_simodelstate_t* state); 315 316 317 initialize 318 ---------- 319 320 This function initializes a given composite model. After this call, the model is in simulation mode. 321 322 PYTHON 323 ^^^^^^ 324 325 Args: 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 331 Returns: 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 338 C 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 346 newSysIdentModel 347 ---------------- 348 349 Creates an empty model for parameter estimation. 350 351 PYTHON 352 ^^^^^^ 353 354 The corresponding Python function is the class constructor. 355 356 Args: 357 :ident: (str) Name of the model instance. 358 359 Returns: 360 :simodel: SysIdent model instance. 361 362 .. code-block:: python 363 364 simodel = OMSysIdent(ident) 365 366 C 367 ^ 368 369 .. code-block:: c 370 371 void* omsi_newSysIdentModel(const char* ident); 372 373 374 setOptions_max_num_iterations 375 ----------------------------- 376 377 Set Ceres solver option `Solver::Options::max_num_iterations`. 378 379 PYTHON 380 ^^^^^^ 381 382 Args: 383 :max_num_iterations: (int) Maximum number of iterations for which the solver should run (default: 25). 384 385 Returns: 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 392 C 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 400 solve 401 ----- 402 403 Solve parameter estimation problem. 404 405 PYTHON 406 ^^^^^^ 407 408 Args: 409 :reporttype: (str) Print report and progress information after call to Ceres solver. 410 Supported report types: `"", "BriefReport", "FullReport"`, where `""` denotes no output. 411 412 Returns: 413 :status: (int) The C-API status code (`oms_status_enu_t`). 414 415 .. code-block:: python 416 417 status = simodel.solve(reporttype) 418 419 C 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.