Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 0.0% 0 / 0 / 532
Functions: 0.0% 0 / 0 / 19
Branches: 0.0% 0 / 0 / 214
OMCompiler/SimulationRuntime/c/optimization/DataManagement/MoveData.c
Line Branch Exec Source
1 /*
2 * This file belongs to the OpenModelica Run-Time System
3 *
4 * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC), c/o Linköpings
5 * universitet, Department of Computer and Information Science, SE-58183 Linköping, Sweden. All rights
6 * reserved.
7 *
8 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE
9 * AGPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. ANY
10 * USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES RECIPIENT'S
11 * ACCEPTANCE OF THE BSD NEW LICENSE OR THE OSMC PUBLIC LICENSE OR THE AGPL
12 * VERSION 3, ACCORDING TO RECIPIENTS CHOICE.
13 *
14 * The OpenModelica software and the OSMC (Open Source Modelica Consortium) Public License
15 * (OSMC-PL) are obtained from OSMC, either from the above address, from the URLs:
16 * http://www.openmodelica.org or https://github.com/OpenModelica/ or
17 * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica distribution. GNU
18 * AGPL version 3 is obtained from: https://www.gnu.org/licenses/licenses.html#GPL. The BSD NEW
19 * License is obtained from: http://www.opensource.org/licenses/BSD-3-Clause.
20 *
21 * This program is distributed WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY
23 * SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF
24 * OSMC-PL.
25 *
26 */
27
28 /*! MoveData.c
29 */
30
31 #include "../../openmodelica_types.h"
32 #include "../../openmodelica.h"
33 #include "../../simulation/arrayIndex.h"
34 #include "../../simulation/options.h"
35 #include "../../simulation/results/simulation_result.h"
36 #include "../../simulation/solver/model_help.h"
37 #include "../../util/real_array.h"
38 #include "../../util/context.h"
39 #include "../../util/omc_file.h"
40 #include "../OptimizerData.h"
41 #include "../OptimizerLocalFunction.h"
42
43 static inline void pickUpDim(OptDataDim * dim, DATA* data, OptDataTime * time);
44 static inline void pickUpTime(OptDataTime * time, OptDataDim * dim, DATA* data, const double preSimTime);
45 static inline void pickUpBounds(OptDataBounds * bounds, OptDataDim * dim, DATA* data);
46 static inline void check_nominal(OptDataBounds * bounds, const double min, const double max,
47 const double nominal, const modelica_boolean set, const int i, const double x0);
48 static inline void calculatedScalingHelper(OptDataBounds * bounds, OptDataTime * time, OptDataDim * dim,OptDataRK * rk);
49
50 static inline void setRKCoeff(OptDataRK *rk, const int np);
51 static inline void printSomeModelInfos(OptDataBounds * bounds, OptDataDim * dim, DATA* data);
52 static inline void pickUpStates(OptData* optdata);
53 static inline void updateDOSystem(OptData * optData, DATA * data, threadData_t *threadData,
54 const int i, const int j, const int index, const int m);
55
56 void setLocalVars(OptData * optData, DATA * data, const double * const vopt, const int i, const int j, const int shift);
57
58 static inline int getNsi(char*, const int, modelica_boolean*);
59 static inline void overwriteTimeGridFile(OptDataTime * time, char* filename, long double c[], const int np, const int nsi);
60 static inline void overwriteTimeGridModel(OptDataTime * time, long double c[], const int np, const int nsi);
61
62 /* pick up model data
63 */
64 ✗ int pickUpModelData(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
65 {
66 ✗ const int nReal = data->modelData->nVariablesReal;
67 ✗ const int nBoolean = data->modelData->nVariablesBoolean;
68 ✗ const int nInteger = data->modelData->nVariablesInteger;
69 ✗ const int nRelations = data->modelData->nRelations;
70
71 int i, j;
72 ✗ OptData *optData = (OptData*) solverInfo->solverData;
73 OptDataDim *dim;
74
75 ✗ pickUpDim(&optData->dim, data, &optData->time);
76 ✗ pickUpBounds(&optData->bounds, &optData->dim, data);
77 ✗ pickUpTime(&optData->time, &optData->dim, data, optData->bounds.preSim);
78 ✗ setRKCoeff(&optData->rk, optData->dim.np);
79 ✗ calculatedScalingHelper(&optData->bounds,&optData->time, &optData->dim, &optData->rk);
80 ✗ messageClose(OMC_LOG_SOLVER); // FIXME what does this belong to?
81
82 dim = &optData->dim;
83
84 ✗ optData->v = (modelica_real***) malloc(dim->nsi*sizeof(modelica_real**));
85 ✗ for(i = 0; i< dim->nsi; ++i){
86 ✗ optData->v[i] = (modelica_real**)malloc(dim->np*sizeof(modelica_real*));
87 ✗ for(j = 0; j<dim->np;++j)
88 ✗ optData->v[i][j] = (modelica_real*)malloc(nReal*sizeof(modelica_real));
89 }
90 ✗ optData->data = data;
91 ✗ optData->threadData = threadData;
92
93 ✗ optData->v0 = (modelica_real*)malloc(nReal*sizeof(modelica_real));
94 ✗ memcpy(optData->v0, data->localData[0]->realVars, nReal*sizeof(modelica_real));
95
96 ✗ pickUpStates(optData);
97
98 ✗ optData->sv0 = (modelica_real*)malloc(dim->nx*sizeof(modelica_real));
99 ✗ for(i = 0; i<dim->nx; ++i)
100 ✗ optData->sv0[i] = optData->v0[i] * optData->bounds.scalF[i];
101
102 ✗ optData->i0 = (modelica_integer*)malloc(nInteger*sizeof(modelica_integer));
103 ✗ memcpy(optData->i0, data->localData[0]->integerVars, nInteger*sizeof(modelica_integer));
104
105 ✗ optData->b0 = (modelica_boolean*)malloc(nBoolean*sizeof(modelica_boolean));
106 ✗ memcpy(optData->b0, data->localData[0]->booleanVars, nBoolean*sizeof(modelica_boolean));
107
108 ✗ optData->re = (modelica_boolean*)malloc(nRelations*sizeof(modelica_boolean));
109 ✗ memcpy(optData->re, data->simulationInfo->relations, nRelations*sizeof(modelica_boolean));
110
111 ✗ optData->i0Pre = (modelica_integer*)malloc(nInteger*sizeof(modelica_integer));
112 ✗ memcpy(optData->i0Pre, data->simulationInfo->integerVarsPre, nInteger*sizeof(modelica_integer));
113
114 ✗ optData->b0Pre = (modelica_boolean*)malloc(nBoolean*sizeof(modelica_boolean));
115 ✗ memcpy(optData->b0Pre, data->simulationInfo->booleanVarsPre, nBoolean*sizeof(modelica_boolean));
116
117 ✗ optData->v0Pre = (modelica_real*)malloc(nReal*sizeof(modelica_real));
118 ✗ memcpy(optData->v0Pre, data->simulationInfo->realVarsPre, nReal*sizeof(modelica_real));
119
120 ✗ optData->rePre = (modelica_boolean*)malloc(nRelations*sizeof(modelica_boolean));
121 ✗ memcpy(optData->rePre, data->simulationInfo->relationsPre, nRelations*sizeof(modelica_boolean));
122
123 ✗ optData->storeR = (modelica_boolean*)malloc(nRelations*sizeof(modelica_boolean));
124 ✗ memcpy(optData->storeR, data->simulationInfo->storedRelations, nRelations*sizeof(modelica_boolean));
125
126 ✗ printSomeModelInfos(&optData->bounds, &optData->dim, data);
127
128 ✗ return 0;
129 }
130
131 /* pick up information(nStates...) from model data to optimizer struct
132 */
133 ✗ static inline void pickUpDim(OptDataDim * dim, DATA* data, OptDataTime * time){
134 char * cflags = NULL;
135 ✗ cflags = (char*)omc_flagValue[FLAG_OPTIMIZER_NP];
136 ✗ if (cflags) {
137 ✗ dim->np = atoi(cflags);
138 ✗ if (dim->np != 1 && dim->np!=3) {
139 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "FLAG_OPTIZER_NP is %i. Currently optimizer support only 1 and 3.\nFLAG_OPTIZER_NP set of 3", dim->np);
140 ✗ dim->np = 3;
141 }
142 } else {
143 ✗ dim->np = 3; /*ToDo*/
144 }
145 ✗ dim->nx = data->modelData->nStates;
146 ✗ dim->nu = data->modelData->nInputVars;
147 ✗ dim->nv = dim->nx + dim->nu;
148 ✗ dim->nc = data->modelData->nOptimizeConstraints;
149 ✗ dim->ncf = data->modelData->nOptimizeFinalConstraints;
150 ✗ dim->nJ = dim->nx + dim->nc;
151 ✗ dim->nJ2 = dim->nJ + 2;
152 ✗ dim->nReal = data->modelData->nVariablesReal;
153
154 ✗ cflags = (char*)omc_flagValue[FLAG_OPTIMIZER_TGRID];
155 ✗ dim->nsi = -1; /* Initialize the data just in case */
156 {
157 /* The model names its time grid by parameter index; the values are read here. */
158 ✗ modelica_integer *tgrid = NULL;
159 modelica_integer i;
160 ✗ data->callback->getTimeGrid(data, &dim->nsi, &tgrid); /* TODO: dim->nsi is long*, expected is int* */
161 ✗ if (dim->nsi > 0) {
162 ✗ time->tt = (modelica_real*) malloc((dim->nsi+1)*sizeof(modelica_real));
163 ✗ for (i = 0; i < dim->nsi+1; ++i) {
164 ✗ time->tt[i] = data->simulationInfo->realParameter[tgrid[i]];
165 }
166 }
167 ✗ free(tgrid);
168 }
169 ✗ time->model_grid = (modelica_boolean)(dim->nsi > 0);
170
171 ✗ if (!time->model_grid) {
172 ✗ dim->nsi = data->simulationInfo->numSteps;
173 }
174
175 ✗ if (cflags) {
176 ✗ dim->nsi = getNsi(cflags, dim->nsi, &dim->exTimeGrid);
177 }
178
179 ✗ dim->nt = dim->nsi*dim->np;
180 ✗ dim->NV = dim->nt*dim->nv;
181 ✗ dim->NRes = dim->nt*dim->nJ + dim->ncf;
182 ✗ dim->index_con = dim->nReal - (dim->nc + dim->ncf);
183 ✗ dim->index_conf = dim->index_con + dim->nc;
184 ✗ assert(dim->nt > 0);
185 ✗ }
186
187
188
189 /* pick up information(startTime, stopTime, dt) from model data to optimizer struct
190 */
191 ✗ static inline void pickUpTime(OptDataTime * time, OptDataDim * dim, DATA* data, const double preSimTime){
192 ✗ const int nsi = dim->nsi;
193 ✗ const int np = dim->np;
194 ✗ const int np1 = np - 1;
195 ✗ long double *c = (long double*)malloc(np * sizeof(long double));
196 ✗ long double *dc = (long double*)malloc(np * sizeof(long double));
197 int i, k;
198 double t;
199 char * cflags = NULL;
200
201 ✗ time->t0 = (long double)fmax(data->simulationInfo->startTime, preSimTime);
202 ✗ time->tf = (long double)data->simulationInfo->stopTime;
203
204 ✗ time->dt = (long double*) malloc((nsi+1)*sizeof(long double));
205 ✗ time->dt[0] = (time->tf - time->t0)/nsi;
206
207 ✗ time->t = (long double**)malloc(nsi*sizeof(long double*));
208 ✗ for(i = 0; i<nsi; ++i)
209 ✗ time->t[i] = (long double*)malloc(np*sizeof(long double));
210 ✗ if(nsi < 1){
211 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Not support numberOfIntervals = %i < 1", nsi);
212 ✗ assert(0);
213 }
214
215 ✗ if(np == 1){
216 ✗ c[0] = 1.0;
217 ✗ }else if(np == 3){
218 ✗ c[0] = 0.15505102572168219018027159252941086080340525193433;
219 ✗ c[1] = 0.64494897427831780981972840747058913919659474806567;
220 ✗ c[2] = 1.00000;
221 }else{
222 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Not support np = %i", np);
223 ✗ assert(0);
224 }
225
226 ✗ for(k = 0; k < np; ++k){
227 ✗ dc[k] = c[k]*time->dt[0];
228 ✗ time->t[0][k] = time->t0 + dc[k];
229 }
230
231 ✗ for(i = 1; i < nsi; ++i){
232 ✗ time->dt[i] = time->dt[i-1];
233 ✗ for(k = 0; k < np; ++k)
234 ✗ time->t[i][k] = time->t[i-1][np1] + dc[k];
235 }
236 ✗ time->t[nsi-1][np1] = time->tf;
237
238 ✗ if(nsi > 1){
239 ✗ i = nsi - 1;
240 ✗ time->dt[nsi-1] = time->t[i][np1] - time->t[i-1][np1];
241 ✗ for(k = 0; k < np; ++k)
242 ✗ time->t[i][k] = time->t[i-1][np1] + c[k]*time->dt[nsi-1];
243 }else
244 ✗ time->dt[1] = time->dt[0];
245
246 ✗ cflags = (char*)omc_flagValue[FLAG_OPTIMIZER_TGRID];
247
248 ✗ if(cflags)
249 ✗ overwriteTimeGridFile(time, cflags, c, np, nsi);
250 ✗ if(time->model_grid)
251 ✗ overwriteTimeGridModel(time, c, np, nsi);
252
253 ✗ free(c);
254 ✗ free(dc);
255 ✗ }
256
257 ✗ static int getNsi(char*filename, const int nsi, modelica_boolean * exTimeGrid){
258 int n = 0, c;
259 FILE * pFile = NULL;
260
261 ✗ *exTimeGrid = 0;
262 ✗ pFile = omc_fopen(filename,"r");
263 ✗ if(pFile == NULL){
264 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "OMC can't find the file %s.", filename);
265 ✗ return nsi;
266 }
267 while(1){
268 ✗ c = fgetc(pFile);
269 ✗ if (c==EOF) break;
270 ✗ if (c=='\n') ++n;
271 }
272 // check if csv file is empty!
273 ✗ if (n == 0){
274 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "time grid file: %s is empty", filename);
275 ✗ fclose(pFile);
276 ✗ return nsi;
277 }
278 ✗ *exTimeGrid = 1;
279 ✗ return n-1;
280 }
281
282 ✗ static inline void overwriteTimeGridFile(OptDataTime * time, char* filename, long double c[], const int np, const int nsi){
283 int i,k;
284 ✗ long double *dc = (long double*)malloc(np * sizeof(long double));
285 ✗ const int np1 = np - 1;
286 double t;
287 FILE * pFile = NULL;
288 ✗ pFile = omc_fopen(filename,"r");
289
290 ✗ fscanf(pFile, "%lf", &t);
291 ✗ time->t0 = t;
292 ✗ fscanf(pFile, "%lf", &t);
293 ✗ time->t[0][np1] = t;
294 ✗ time->dt[0] = time->t[0][np1] - time->t0;
295
296 ✗ if(time->dt[0] <= 0){
297 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "read time grid from file fail!");
298 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "line %i: %g <= %g",0, (double)time->t[0][np1], (double)time->t0);
299 ✗ EXIT(0);
300 }
301
302
303 ✗ for(k = 0; k < np; ++k){
304 ✗ dc[k] = c[k]*time->dt[0];
305 ✗ time->t[0][k] = time->t0 + dc[k];
306 }
307
308 ✗ for(i=1;i<nsi;++i){
309 ✗ fscanf(pFile, "%lf", &t);
310 ✗ time->t[i][np1] = t;
311 ✗ time->dt[i] = time->t[i][np1] - time->t[i-1][np1];
312
313 ✗ for(k = 0; k < np; ++k){
314 ✗ dc[k] = c[k]*time->dt[i];
315 ✗ time->t[i][k] = time->t[i-1][np1] + dc[k];
316 }
317
318 ✗ if(time->dt[i] <= 0){
319 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "read time grid");
320 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "line %i/%i: %g <= %g",i, nsi, (double)time->t[i][np1], (double)time->t[i-1][np1]);
321 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "failed!");
322 ✗ EXIT(0);
323 }
324
325 }
326 ✗ time->tf = time->t[nsi-1][np1];
327 ✗ fclose(pFile);
328 ✗ free(dc);
329 ✗ }
330
331 ✗ int cmp_modelica_real(const void *v1, const void *v2) {
332 ✗ return (*(modelica_real*)v1 - *(modelica_real*)v2);
333 }
334
335 ✗ static inline void overwriteTimeGridModel(OptDataTime * time, long double c[], const int np, const int nsi){
336 int i,k;
337
338 ✗ time->t0 = time->tt[0];
339 ✗ time->tf = time->tt[nsi];
340
341 ✗ qsort((void*) time->tt, nsi+1, sizeof(modelica_real), &cmp_modelica_real);
342
343 ✗ for(i = 0; i<nsi; ++i){
344 ✗ time->dt[i] = time->tt[i+1] - time->tt[i];
345 ✗ for(k=0; k<np; ++k){
346 ✗ time->t[i][k] = time->tt[i] + c[k]*time->dt[i];
347 /*printf("\nt[%i][%i] = %g",i,k,(double)time->t[i][k]);*/
348 }
349 }
350
351 ✗ free(time->tt);
352 ✗ }
353
354 /* pick up information(startTime, stopTime, dt) from model data to optimizer struct
355 */
356 ✗ static inline void pickUpBounds(OptDataBounds * bounds, OptDataDim * dim, DATA* data){
357 char ** inputName;
358 double min, max, nominal, x0;
359 double *umin, *umax, *unom;
360 modelica_boolean nominalWasSet;
361 modelica_boolean * nominalWasSetInput;
362
363 ✗ const int nx = dim->nx;
364 ✗ const int nv = dim->nv;
365 ✗ const int nu = dim->nu;
366 ✗ const int nt = dim->nt;
367 ✗ const int NV = dim->NV;
368
369 long double tmp;
370
371 int i, j;
372
373 ✗ dim->inputName = (char**) malloc(nv*sizeof(char*));
374 ✗ bounds->vnom = malloc(nv*sizeof(double));
375 ✗ bounds->scalF = malloc(nv*sizeof(long double));
376
377 ✗ bounds->vmin = malloc(nv*sizeof(double));
378 ✗ bounds->vmax = malloc(nv*sizeof(double));
379
380 ✗ bounds->u0 = malloc(nu*sizeof(double));
381
382 ✗ nominalWasSetInput = (modelica_boolean*)malloc(nv*sizeof(modelica_boolean));
383 inputName = dim->inputName;
384
385 ✗ umin = bounds->vmin + nx;
386 ✗ umax = bounds->vmax + nx;
387 ✗ unom = bounds->vnom + nx;
388
389 ✗ data->callback->pickUpBoundsForInputsInOptimization(data,umin, umax, unom, nominalWasSetInput, inputName, bounds->u0, &bounds->preSim);
390
391 ✗ for(i = 0; i < nx; ++i){
392 ✗ min = getMinFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, i);
393 ✗ max = getMaxFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, i);
394 ✗ nominal = getNominalFromScalarIdx(data->simulationInfo, data->modelData, VAR_KIND_VARIABLE, i);
395 ✗ nominalWasSet = data->modelData->realVarsData[i].attribute.useNominal;
396 ✗ x0 = data->localData[1]->realVars[i];
397
398 ✗ check_nominal(bounds, min, max, nominal, nominalWasSet, i, x0);
399 ✗ array_index_t* revIndex = &data->simulationInfo->realVarsReverseIndex[i];
400 ✗ put_real_element(bounds->vnom[i], revIndex->dim_idx, &data->modelData->realVarsData[revIndex->array_idx].attribute.nominal);
401 ✗ bounds->scalF[i] = 1.0/bounds->vnom[i];
402 ✗ bounds->vmin[i] = min * bounds->scalF[i];
403 ✗ bounds->vmax[i] = max * bounds->scalF[i];
404
405 }
406 ✗ for(j=0; i<dim->nv; ++i,++j){
407
408 ✗ bounds->u0[j] = fmin(fmax(bounds->u0[j], umin[j]), umax[j]);
409 ✗ check_nominal(bounds, umin[j], umax[j], unom[j], nominalWasSetInput[j], i, fabs(bounds->u0[j]));
410
411 ✗ bounds->scalF[i] = 1.0 / bounds->vnom[i];
412 ✗ bounds->vmin[i] *= bounds->scalF[i];
413 ✗ bounds->vmax[i] *= bounds->scalF[i];
414
415 }
416 ✗ free(nominalWasSetInput);
417
418 ✗ bounds->Vmin = malloc(NV*sizeof(double));
419 ✗ bounds->Vmax = malloc(NV*sizeof(double));
420
421 ✗ for(i = 0, j = 0; i < nt ; ++i, j += nv){
422 ✗ memcpy(bounds->Vmin + j, bounds->vmin, nv*sizeof(double));
423 ✗ memcpy(bounds->Vmax + j, bounds->vmax, nv*sizeof(double));
424 }
425
426 ✗ }
427
428 /*!
429 * heuristic for nominal value
430 * author: Vitalij Ruge
431 **/
432 ✗ static inline void check_nominal(OptDataBounds * bounds, const double min, const double max,
433 const double nominal, const modelica_boolean set, const int i, const double x0){
434
435 ✗ if(set){
436 ✗ bounds->vnom[i] = fmax(fabs(nominal),1e-16);
437 }else{
438 double amax, amin;
439
440 ✗ amax = fabs(max);
441 ✗ amin = fabs(min);
442
443 ✗ bounds->vnom[i] = fmax(amax,amin);
444
445 ✗ if(bounds->vnom[i] > 1e12){
446 ✗ double tmp = fmin(amax,amin);
447 ✗ double ax0 = fabs(x0);
448 ✗ bounds->vnom[i] = (tmp < 1e12) ? fmax(tmp,ax0) : 1.0 + ax0;
449 }
450
451 ✗ bounds->vnom[i] = fmax(bounds->vnom[i], 1e-16);
452 }
453 ✗ }
454
455 /*!
456 * calculated helper vars for scaling
457 * author: Vitalij Ruge
458 **/
459 ✗ static inline void calculatedScalingHelper(OptDataBounds * bounds, OptDataTime * time, OptDataDim * dim, OptDataRK *rk){
460 ✗ const int nx = dim->nx;
461 ✗ const int nsi = dim->nsi;
462 ✗ const int np = dim->np;
463
464 int i, j, k, l;
465 ✗ assert(nsi > 0);
466 ✗ bounds->scaldt = (long double**)malloc(nsi*sizeof(long double*));
467 ✗ for(i = 0; i < nsi; ++i)
468 ✗ bounds->scaldt[i] = (long double*) malloc(nx*sizeof(long double));
469
470 ✗ for(i = 0; i < nsi; ++i)
471 ✗ for(j = 0; j < nx; ++j){
472 ✗ bounds->scaldt[i][j] = bounds->scalF[j]*time->dt[i];
473 }
474
475 ✗ bounds->scalb = (long double**)malloc(nsi*sizeof(long double*));
476 ✗ for(i = 0; i < nsi; ++i){
477 ✗ bounds->scalb[i] = (long double*)malloc(np*sizeof(long double));
478 ✗ for(j = 0; j < np; ++j){
479 ✗ bounds->scalb[i][j] = time->dt[i]*rk->b[j];
480 }
481 }
482 ✗ }
483
484 /*!
485 * set RK coeffs
486 * author: Vitalij Ruge
487 **/
488 static inline void setRKCoeff(OptDataRK *rk, const int np){
489
490 ✗ if(np == 3){
491
492 ✗ rk->a[0][0] = 4.1393876913398137178367408896470696703591369767880;
493 ✗ rk->a[0][1] = 3.2247448713915890490986420373529456959829737403284;
494 ✗ rk->a[0][2] = 1.1678400846904054949240412722156950122337492313015;
495 ✗ rk->a[0][3] = 0.25319726474218082618594241992157103785758599484179;
496
497 ✗ rk->a[1][0] = 1.7393876913398137178367408896470696703591369767880;
498 ✗ rk->a[1][1] = 3.5678400846904054949240412722156950122337492313015;
499 ✗ rk->a[1][2] = 0.7752551286084109509013579626470543040170262596716;
500 ✗ rk->a[1][3] = 1.0531972647421808261859424199215710378575859948418;
501
502 ✗ rk->a[2][0] = 3.0;
503 ✗ rk->a[2][1] = 5.5319726474218082618594241992157103785758599484179;
504 ✗ rk->a[2][2] = 7.5319726474218082618594241992157103785758599484179;
505 ✗ rk->a[2][3] = 5.0;
506
507 ✗ rk->b[0] = 0.37640306270046727505007544236928079466761256998175;
508 ✗ rk->b[1] = 0.51248582618842161383881344651960809422127631890713;
509 ✗ rk->b[2] = 1 - (rk->b[0] + rk->b[1]);
510
511 ✗ }else if(np == 1){
512 ✗ rk->a[0][0] = 1.000;
513 ✗ rk->b[0] = rk->a[0][0];
514 }
515 }
516
517 /*!
518 * print some model infos
519 * author: Vitalij Ruge
520 **/
521 ✗ static inline void printSomeModelInfos(OptDataBounds * bounds, OptDataDim * dim, DATA* data)
522 {
523 ✗ const int nx = dim->nx;
524 ✗ const int nc = dim->nc;
525 ✗ const int nv = dim->nv;
526
527 double *umin, *umax, *unom, *u0;
528 double *xmin, *xmax, *xnom;
529
530 char buffer[200];
531
532 char ** inputName;
533 int i,j,k;
534
535 ✗ inputName = dim->inputName;
536
537 ✗ umin = bounds->vmin + nx;
538 ✗ umax = bounds->vmax + nx;
539 ✗ unom = bounds->vnom + nx;
540 ✗ u0 = bounds->u0;
541
542 xmin = bounds->vmin;
543 xmax = bounds->vmax;
544 xnom = bounds->vnom;
545
546 printf("\nOptimizer Variables");
547 printf("\n========================================================");
548
549 ✗ for(i = 0; i < nx; ++i){
550
551 ✗ if (xmin[i] > -1e20) {
552 ✗ sprintf(buffer, ", min = %g", real_get(data->modelData->realVarsData[i].attribute.min, 0));
553 }
554 else {
555 sprintf(buffer, ", min = -Inf");
556 }
557
558 ✗ printf("\nState[%i]:%s(start = %g, nominal = %g%s",
559 i,
560 ✗ data->modelData->realVarsData[i].info.name,
561 ✗ real_get(data->modelData->realVarsData[i].attribute.start, 0),
562 ✗ xnom[i],
563 buffer);
564
565 ✗ if(xmax[i] < 1e20)
566 ✗ sprintf(buffer, ", max = %g", real_get(data->modelData->realVarsData[i].attribute.max, 0));
567 else
568 sprintf(buffer, ", max = +Inf");
569
570 printf("%s",buffer);
571 ✗ printf(", init = %g)", data->localData[1]->realVars[i]);
572 }
573
574 ✗ for(k = 0; i < nv; ++i, ++k){
575
576 ✗ if (umin[k] > -1e20)
577 ✗ sprintf(buffer, ", min = %g", umin[k]*unom[k]);
578 else
579 sprintf(buffer, ", min = -Inf");
580
581 ✗ printf("\nInput[%i]:%s(start = %g, nominal = %g%s",i, inputName[k], u0[k], unom[k], buffer);
582
583 ✗ if(umax[k] < 1e20)
584 ✗ sprintf(buffer, ", max = %g", umax[k]*unom[k]);
585 else
586 sprintf(buffer, ", max = +Inf");
587
588 printf("%s)",buffer);
589 }
590 printf("\n--------------------------------------------------------");
591 printf("\nnumber of nonlinear constraints: %i", nc);
592 printf("\n========================================================\n");
593
594 ✗ }
595
596
597 /*!
598 * write results in result file
599 * author: Vitalij Ruge
600 **/
601 ✗ void res2file(OptData *optData, SOLVER_INFO* solverInfo, double *vopt){
602 ✗ const int nu = optData->dim.nu;
603 ✗ const int nx = optData->dim.nx;
604 ✗ const int nv = optData->dim.nv;
605 ✗ const int nsi = optData->dim.nsi;
606 ✗ const int np = optData->dim.np;
607 ✗ const int nReal = optData->dim.nReal;
608 ✗ const int nBoolean = optData->data->modelData->nVariablesBoolean;
609 const int nInteger = optData->data->modelData->nVariablesInteger;
610 const int nRelations = optData->data->modelData->nRelations;
611 ✗ const int nvnp = nv*np;
612 ✗ long double *a = (long double*)malloc(np * sizeof(long double));
613 ✗ modelica_real *** v = optData->v;
614 float tmp_u;
615
616 int i,j,k, ii, jj;
617 char buffer[4096];
618 DATA * data = optData->data;
619 ✗ threadData_t *threadData = optData->threadData;
620 ✗ SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
621
622 ✗ FILE * pFile = optData->pFile;
623 ✗ double *vnom = optData->bounds.vnom;
624 ✗ long double **t = optData->time.t;
625 ✗ long double t0 = optData->time.t0;
626 long double tmpv;
627
628 ✗ if(np == 3){
629 ✗ a[0] = 1.5580782047249223824319753706862790293163070736617;
630 ✗ a[1] = -0.89141153805825571576530870401961236264964040699507;
631 ✗ a[2] = 0.33333333333333333333333333333333333333333333333333;
632 ✗ }else if(np == 1){
633 ✗ a[0] = 1.000;
634 }else{
635 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Not support np = %i", np);
636 ✗ assert(0);
637 }
638
639 ✗ optData2ModelData(optData, vopt, 0);
640
641 /******************/
642 ✗ fprintf(pFile, "%lf ",(double)t0);
643
644 ✗ for(i=0,j = nx; i < nu; ++i,++j){
645 ✗ for(k = 0, tmpv = 0.0; k < np; ++k){
646 ✗ tmpv += a[k]*vopt[k*nv + j];
647 }
648 ✗ tmpv = fmin(fmax(tmpv,optData->bounds.vmin[j]),optData->bounds.vmax[j]);
649 ✗ data->simulationInfo->inputVars[i] = (double)tmpv*vnom[j];
650 ✗ fprintf(pFile, "%lf ", (float)data->simulationInfo->inputVars[i]);
651 }
652 fprintf(pFile, "%s", "\n");
653 /******************/
654 ✗ copy_initial_values(optData, data);
655 /******************/
656 ✗ solverInfo->currentTime = (double)t0;
657 ✗ sData->timeValue = solverInfo->currentTime;
658
659 /*updateDiscreteSystem(data);*/
660 ✗ data->callback->input_function(data, threadData);
661 /*data->callback->functionDAE(data);*/
662 ✗ updateDiscreteSystem(data, threadData);
663
664 ✗ sim_result.emit(&sim_result, data, threadData);
665 /******************/
666
667 ✗ for(ii = 0; ii < nsi; ++ii){
668 ✗ for(jj = 0; jj < np; ++jj){
669 /******************/
670 ✗ memcpy(sData->realVars, v[ii][jj], nReal*sizeof(modelica_real));
671 /******************/
672 ✗ fprintf(pFile, "%lf ",(double)t[ii][jj]);
673 ✗ for(i = 0; i < nu; ++i){
674 ✗ tmp_u = (float)(vopt[ii*nvnp+jj*nv+nx+i]*vnom[i + nx]);
675 ✗ fprintf(pFile, "%lf ", tmp_u);
676 }
677 fprintf(pFile, "%s", "\n");
678 /******************/
679 ✗ solverInfo->currentTime = (double)t[ii][jj];
680 ✗ sData->timeValue = solverInfo->currentTime;
681 ✗ sim_result.emit(&sim_result, data, threadData);
682 }
683 }
684 ✗ fclose(pFile);
685 ✗ free(a);
686 ✗ }
687
688
689 ✗ void copy_initial_values(OptData * optData, DATA* data){
690 ✗ const int nBoolean = optData->data->modelData->nVariablesBoolean ;
691 ✗ const int nInteger = optData->data->modelData->nVariablesInteger;
692 ✗ const int nReal = optData->dim.nReal;
693 ✗ const int nRelations = optData->data->modelData->nRelations;
694
695 ✗ memcpy(data->localData[0]->realVars, optData->v0, nReal*sizeof(modelica_real));
696 ✗ memcpy(data->localData[0]->integerVars, optData->i0, nInteger*sizeof(modelica_integer));
697 ✗ memcpy(data->localData[0]->booleanVars, optData->b0, nBoolean*sizeof(modelica_boolean));
698 ✗ memcpy(data->simulationInfo->integerVarsPre, optData->i0Pre, nInteger*sizeof(modelica_integer));
699 ✗ memcpy(data->simulationInfo->booleanVarsPre, optData->b0Pre, nBoolean*sizeof(modelica_boolean));
700 ✗ memcpy(data->simulationInfo->realVarsPre, optData->v0Pre, nReal*sizeof(modelica_real));
701 ✗ memcpy(data->simulationInfo->relationsPre, optData->rePre, nRelations*sizeof(modelica_boolean));
702 ✗ memcpy(data->simulationInfo->relations, optData->re, nRelations*sizeof(modelica_boolean));
703 ✗ memcpy(data->simulationInfo->storedRelations, optData->storeR, nRelations*sizeof(modelica_boolean));
704
705 ✗ }
706
707 /*!
708 * transfer optimizer data to model data
709 * author: Vitalij Ruge
710 **/
711 ✗ void optData2ModelData(OptData *optData, double *vopt, const int index){
712 ✗ const int nv = optData->dim.nv;
713 ✗ const int nsi = optData->dim.nsi;
714 ✗ const int np = optData->dim.np;
715
716 modelica_real * realVars[3];
717 ✗ modelica_real * tmpVars[2] = {NULL, NULL};
718
719 int i, j, k, shift, l;
720 ✗ DATA * data = optData->data;
721 ✗ const int * indexBC = optData->s.indexABCD + 3;
722 ✗ threadData_t *threadData = optData->threadData;
723
724 ✗ for(l = 0; l < 3; ++l)
725 ✗ realVars[l] = data->localData[l]->realVars;
726
727 ✗ for(l = 0; l< 2; ++l){
728 ✗ if(optData->s.matrix[l])
729 ✗ tmpVars[l] = data->simulationInfo->analyticJacobians[indexBC[l]].tmpVars;
730 }
731 ✗ copy_initial_values(optData, data);
732
733 ✗ for(i = 0, shift = 0; i < nsi-1; ++i){
734 ✗ for(j = 0; j < np; ++j, shift += nv){
735 ✗ setLocalVars(optData, data, vopt, i, j, shift);
736 ✗ updateDOSystem(optData, data, threadData, i, j, index, 2);
737 }
738 }
739
740 ✗ for(j = 0; j < np-1; ++j, shift += nv){
741 ✗ setLocalVars(optData, data, vopt, i, j, shift);
742 ✗ updateDOSystem(optData, data, threadData, i, j, index, 2);
743 }
744 ✗ setLocalVars(optData, data, vopt, i, j, shift);
745 ✗ updateDOSystem(optData, data, threadData, i, j, index, 3);
746
747 /*terminal constraint(s)*/
748 ✗ if(index){
749 ✗ if(optData->s.matrix[3])
750 ✗ diffSynColoredOptimizerSystemF(optData, optData->Jf);
751 }
752
753 ✗ for(l = 0; l < 3; ++l)
754 ✗ data->localData[l]->realVars = realVars[l];
755
756 ✗ for(l = 0; l< 2; ++l)
757 ✗ if(optData->s.matrix[l])
758 ✗ data->simulationInfo->analyticJacobians[indexBC[l]].tmpVars = tmpVars[l];
759
760 ✗ }
761
762
763 /*!
764 * helper optData2ModelData
765 * author: Vitalij Ruge
766 **/
767 ✗ static inline void updateDOSystem(OptData * optData, DATA * data, threadData_t *threadData,
768 const int i, const int j, const int index, const int m){
769
770 /* try */
771 ✗ optData->scc = 0;
772 #if !defined(OMC_EMCC)
773 ✗ OMC_TRY_INTERNAL(simulationJumpBuffer)
774 #endif
775 ✗ data->callback->input_function(data, optData->threadData);
776 ✗ updateDiscreteSystem(data, optData->threadData);
777
778 ✗ if(index){
779 ✗ diffSynColoredOptimizerSystem(optData, optData->J[i][j], i, j, m);
780 }
781 ✗ optData->scc = 1;
782 #if !defined(OMC_EMCC)
783 ✗ OMC_CATCH_INTERNAL(simulationJumpBuffer)
784 #endif
785 ✗ }
786
787 /*!
788 * helper optData2ModelData
789 * author: Vitalij Ruge
790 **/
791 ✗ void setLocalVars(OptData * optData, DATA * data, const double * const vopt,
792 const int i, const int j, const int shift){
793 short l;
794 int k;
795
796 ✗ const int * indexBC = optData->s.indexABCD + 3;
797 OptDataDim * dim = &optData->dim;
798 ✗ const modelica_real * vnom = optData->bounds.vnom;
799 ✗ const int nx = optData->dim.nx;
800 ✗ const int nv = optData->dim.nv;
801
802 /* try to init discrete real variables with pre value */
803 ✗ memcpy(optData->v[i][j], data->simulationInfo->realVarsPre, optData->dim.nReal*sizeof(modelica_real));
804 ✗ for(l = 0; l < 3; ++l){
805 ✗ data->localData[l]->realVars = optData->v[i][j];
806 ✗ data->localData[l]->timeValue = (modelica_real) optData->time.t[i][j];
807 }
808
809 ✗ for(l = 0; l < 2; ++l)
810 ✗ if(optData->s.matrix[l])
811 ✗ data->simulationInfo->analyticJacobians[indexBC[l]].tmpVars = dim->analyticJacobians_tmpVars[l][i][j];
812
813 ✗ for(k = 0; k < nx; ++k)
814 ✗ data->localData[0]->realVars[k] = vopt[shift + k]*vnom[k];
815
816 ✗ for(; k <nv; ++k){
817 ✗ data->simulationInfo->inputVars[k-nx] = (modelica_real) vopt[shift + k]*vnom[k];
818 }
819
820 ✗ }
821
822
823 /*
824 * function calculates a symbolic colored jacobian matrix of the optimization system
825 * authors: Willi Braun, Vitalij Ruge
826 */
827 ✗ void diffSynColoredOptimizerSystem(OptData *optData, modelica_real **J, const int m, const int n, const int index){
828 ✗ DATA * data = optData->data;
829 ✗ threadData_t *threadData = optData->threadData;
830 int i,j,l,ii, ll;
831
832 ✗ const int h_index = optData->s.indexABCD[index];
833 ✗ JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[h_index]);
834 ✗ const long double * scaldt = optData->bounds.scaldt[m];
835 ✗ const unsigned int * const cC = jacobian->sparsePattern->colorCols;
836 ✗ const unsigned int * const lindex = jacobian->sparsePattern->leadindex;
837 ✗ const int nx = jacobian->sizeCols;
838 ✗ const int Cmax = jacobian->sparsePattern->maxColors + 1;
839 ✗ const int dnx = optData->dim.nx;
840 ✗ const int dnxnc = optData->dim.nJ;
841 ✗ const modelica_real * const resultVars = jacobian->resultVars;
842 ✗ const unsigned int * const sPindex = jacobian->sparsePattern->index;
843 ✗ long double scalb = optData->bounds.scalb[m][n];
844
845 ✗ const int * index_J = (index == 3)? optData->s.indexJ3 : optData->s.indexJ2;
846 ✗ const int nJ1 = optData->dim.nJ + 1;
847
848 ✗ modelica_real **sV = optData->s.seedVec[index];
849 /* The optimizer lends the Jacobian a seed vector of its own per colour. The
850 Jacobian owns seedVars and frees it, so give its own back. */
851 ✗ modelica_real * const ownSeedVars = jacobian->seedVars;
852
853 /* set symbolic jacobian context to reuse the matrix and the factorization in every column */
854 ✗ setContext(data, data->localData[0]->timeValue, CONTEXT_SYM_JACOBIAN);
855
856 ✗ if (jacobian->constantEqns != NULL) {
857 ✗ jacobian->constantEqns(data, threadData, jacobian, NULL);
858 }
859
860 ✗ for(i = 1; i < Cmax; ++i){
861 ✗ jacobian->seedVars = sV[i];
862
863 ✗ if(index == 2){
864 ✗ data->callback->functionJacB_column(data, threadData, jacobian, NULL);
865 ✗ }else if(index == 3){
866 ✗ data->callback->functionJacC_column(data, threadData, jacobian, NULL);
867 }else
868 ✗ assert(0);
869
870 ✗ increaseJacContext(data);
871
872 ✗ for(ii = 0; ii < nx; ++ii){
873 ✗ if(cC[ii] == i){
874 ✗ for(j = lindex[ii]; j < lindex[ii + 1]; ++j){
875 ✗ ll = sPindex[j];
876 ✗ l = index_J[ll];
877 ✗ if(l < dnx){
878 ✗ J[l][ii] = (modelica_real) resultVars[ll] * scaldt[l];
879 ✗ }else if(l < dnxnc){
880 ✗ J[l][ii] = (modelica_real) resultVars[ll];
881 ✗ }else if(l == optData->dim.nJ && optData->s.lagrange){
882 ✗ J[l][ii] = (modelica_real) resultVars[ll]* scalb;
883 ✗ }else if(l == nJ1 && optData->s.mayer){
884 ✗ J[l][ii] = (modelica_real) resultVars[ll];
885 }
886 }
887 }
888
889 }
890 }
891 ✗ jacobian->seedVars = ownSeedVars;
892 /* set context for the start values extrapolation of non-linear algebraic loops */
893 ✗ unsetContext(data);
894 ✗ }
895
896 ✗ void diffSynColoredOptimizerSystemF(OptData *optData, modelica_real **J){
897 ✗ if(optData->dim.ncf > 0){
898 ✗ DATA * data = optData->data;
899 ✗ threadData_t *threadData = optData->threadData;
900 int i,j,l,ii, ll;
901 const int index = 4;
902 ✗ const int h_index = optData->s.indexABCD[index];
903 ✗ JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[h_index]);
904 ✗ const unsigned int * const cC = jacobian->sparsePattern->colorCols;
905 ✗ const unsigned int * const lindex = jacobian->sparsePattern->leadindex;
906 ✗ const int nx = jacobian->sizeCols;
907 ✗ const int Cmax = jacobian->sparsePattern->maxColors + 1;
908 ✗ const modelica_real * const resultVars = jacobian->resultVars;
909 ✗ const unsigned int * const sPindex = jacobian->sparsePattern->index;
910
911 ✗ modelica_real **sV = optData->s.seedVec[index];
912 /* See diffSynColoredOptimizerSystem: seedVars is the Jacobian's to free. */
913 ✗ modelica_real * const ownSeedVars = jacobian->seedVars;
914
915 /* set symbolic jacobian context to reuse the matrix and the factorization in every column */
916 ✗ setContext(data, data->localData[0]->timeValue, CONTEXT_SYM_JACOBIAN);
917
918 ✗ if (jacobian->constantEqns != NULL) {
919 ✗ jacobian->constantEqns(data, threadData, jacobian, NULL);
920 }
921
922 ✗ for(i = 1; i < Cmax; ++i){
923 ✗ jacobian->seedVars = sV[i];
924
925 ✗ data->callback->functionJacD_column(data, threadData, jacobian, NULL);
926
927 ✗ increaseJacContext(data);
928
929 ✗ for(ii = 0; ii < nx; ++ii){
930 ✗ if(cC[ii] == i){
931 ✗ for(j = lindex[ii]; j < lindex[ii + 1]; ++j){
932 ✗ ll = sPindex[j];
933 ✗ J[ll][ii] = resultVars[ll];
934 }
935 }
936 }
937 }
938 ✗ jacobian->seedVars = ownSeedVars;
939 /* set context for the start values extrapolation of non-linear algebraic loops */
940 ✗ unsetContext(data);
941 }
942 ✗ }
943
944 /*!
945 * pick up start values from csv for states
946 * author: Vitalij Ruge
947 **/
948 ✗ static inline void pickUpStates(OptData* optData){
949 char* cflags;
950 ✗ cflags = (char*)omc_flagValue[FLAG_INPUT_FILE_STATES];
951
952 ✗ if(cflags){
953 FILE * pFile = NULL;
954 ✗ pFile = omc_fopen(cflags,"r");
955
956 ✗ if(pFile == NULL){
957 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "OMC can't find the file %s.",cflags);
958 }else{
959 int c, n = 0;
960 modelica_boolean b;
961 while(1){
962 ✗ c = fgetc(pFile);
963 ✗ if (c==EOF) break;
964 ✗ if (c=='\n') ++n;
965 }
966 // check if csv file is empty!
967 ✗ if(n == 0){
968 ✗ fprintf(stderr, "External input file: %s is empty!\n",cflags); fflush(NULL);
969 ✗ EXIT(1);
970 }else{
971 int i, j;
972 double start_value;
973 char buffer[200];
974 ✗ rewind(pFile);
975 ✗ for(i =0; i< n; ++i){
976 ✗ fscanf(pFile, "%199s", buffer);
977 ✗ if (fscanf(pFile, "%lf", &start_value) <= 0) continue;
978
979 ✗ for(j = 0, b = 0; j < optData->dim.nReal; ++j){
980 ✗ if(!strcmp(optData->data->modelData->realVarsData[j].info.name, buffer)){
981 ✗ optData->data->localData[0]->realVars[j] = start_value;
982 ✗ optData->data->localData[1]->realVars[j] = start_value;
983 ✗ optData->data->localData[2]->realVars[j] = start_value;
984 ✗ optData->v0[i] = start_value;
985 b = 1;
986 ✗ continue;
987 }
988 }
989 ✗ if(!b)
990 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "it was impossible to set %s.start %g", buffer,start_value);
991 else
992 ✗ printf("\n[%i]set %s.start %g", i, buffer,start_value);
993
994 }
995 ✗ fclose(pFile);
996 printf("\n");
997 /*update system*/
998 ✗ optData->data->callback->input_function(optData->data, optData->threadData);
999 /*optData->data->callback->functionDAE(optData->data);*/
1000 ✗ updateDiscreteSystem(optData->data, optData->threadData);
1001 }
1002 }
1003 }
1004 ✗ }
1005