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 / 325
Functions: 0.0% 0 / 0 / 10
Branches: 0.0% 0 / 0 / 451
OMCompiler/SimulationRuntime/c/linearization/linearize.cpp
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 #include "util/omc_error.h"
29 #include "util/omc_file.h"
30 #include "simulation_data.h"
31 #include "openmodelica_func.h"
32 #include "simulation/arrayIndex.h"
33 #include "simulation/solver/external_input.h"
34 #include "simulation/options.h"
35 #include "simulation/solver/model_help.h"
36 #include "linearize.h"
37 #include <iostream>
38 #include <sstream>
39 #include <string>
40
41 using namespace std;
42
43 ✗ static string array2string(double* array, int row, int col, DATA *data)
44 {
45 int i=0;
46 int j=0;
47 ✗ ostringstream retVal(ostringstream::out);
48 retVal.precision(16);
49 ✗ for(i=0; i<row; i++)
50 {
51 int k = i;
52 ✗ for(j=0; j<col-1; j++)
53 {
54 ✗ if (data->modelData->linearizationDumpLanguage == OMC_LINEARIZE_DUMP_LANGUAGE_JULIA)
55 {
56 ✗ retVal << array[k] << " "; // Julia matrix accepts space as separators
57 }
58 else
59 {
60 ✗ retVal << array[k] << ", ";
61 }
62 ✗ k += row;
63 }
64 ✗ if(col > 0)
65 {
66 ✗ retVal << array[k];
67 }
68 ✗ if((i+1 != row) && (col != 0))
69 {
70 ✗ retVal << ";\n\t";
71 }
72 }
73 ✗ return retVal.str();
74 ✗ }
75
76 ✗ static string array2PythonString(double* array, int row, int col)
77 {
78 int i=0;
79 int j=0;
80 ✗ ostringstream retVal(ostringstream::out);
81 ✗ if (row == 0 || col == 0)
82 {
83 ✗ retVal << "[]\n";
84 return retVal.str();
85 }
86
87 retVal.precision(16);
88 ✗ retVal << "[";
89 ✗ for(i=0; i<row; i++)
90 {
91 int k = i;
92 ✗ retVal << "[";
93 ✗ for(j=0; j<col-1; j++)
94 {
95 ✗ retVal << array[k] << ", ";
96 ✗ k += row;
97 }
98 ✗ if(col > 0)
99 {
100 ✗ retVal << array[k];
101 }
102 ✗ if((i+1 != row) && (col != 0))
103 {
104 ✗ retVal << "],\n\t";
105 }
106 }
107 ✗ retVal << "]]\n";
108
109 return retVal.str();
110 ✗ }
111
112 extern "C" {
113
114 ✗ int functionODE_residual(DATA* data, threadData_t *threadData, double *dx, double *dy, double *dz)
115 {
116 long i;
117
118 /* debug */
119 /* printCurrentStatesVector(OMC_LOG_JAC, y, data, data->localData[0]->timeValue); */
120
121 /* read input vars */
122 ✗ externalInputUpdate(data);
123 ✗ data->callback->input_function(data, threadData);
124
125 /* eval input vars */
126 ✗ data->callback->functionODE(data, threadData);
127
128 /* eval algebraic vars */
129 ✗ data->callback->functionAlgebraics(data, threadData);
130
131 /* eval output vars */
132 ✗ data->callback->output_function(data, threadData);
133
134 /* get the difference between the temp_xd(=localData->statesDerivatives)
135 and xd(=statesDerivativesBackup) */
136 ✗ for(i=0; i < data->modelData->nStates; i++)
137 {
138 ✗ dx[i] = data->localData[0]->realVars[data->modelData->nStates + i];
139 }
140 ✗ for(i=0; i < data->modelData->nOutputVars; i++)
141 {
142 ✗ dy[i] = data->simulationInfo->outputVars[i];
143 }
144 ✗ if(dz){
145 ✗ for(i=0; i < (data->modelData->nVariablesReal - 2*data->modelData->nStates); i++)
146 {
147 ✗ dz[i] = data->localData[0]->realVars[2*data->modelData->nStates + i];
148 }
149 }
150
151 ✗ return 0;
152 }
153
154 /* Calculate the jacobian matrix by numerical finite difference */
155 ✗ int functionJacAC_num(DATA* data, threadData_t *threadData, double *matrixA, double *matrixC, double *matrixCz)
156 {
157 ✗ const double delta_h = numericalDifferentiationDeltaXlinearize;
158 double delta_hh;
159 double xsave;
160
161 double* x;
162
163 int i,j,k;
164
165 int do_data_recovery = 0;
166
167 ✗ int size_A = data->modelData->nStates;
168 ✗ int size_C = data->modelData->nOutputVars;
169 ✗ int size_z = data->modelData->nVariablesReal - 2*data->modelData->nStates;
170
171 ✗ double* x0 = (double*)calloc(size_A,sizeof(double));
172 ✗ double* y0 = (double*)calloc(size_C,sizeof(double));
173 ✗ double* x1 = (double*)calloc(size_A,sizeof(double));
174 ✗ double* y1 = (double*)calloc(size_C,sizeof(double));
175 double* z0 = 0;
176 double* z1 = 0;
177 ✗ double *xScaling = (double*)calloc(size_A,sizeof(double));
178
179 ✗ assertStreamPrint(threadData,0!=x0,"calloc failed");
180 ✗ assertStreamPrint(threadData,0!=y0,"calloc failed");
181 ✗ assertStreamPrint(threadData,0!=x1,"calloc failed");
182 ✗ assertStreamPrint(threadData,0!=y1,"calloc failed");
183
184 ✗ if(matrixCz){
185 do_data_recovery = 1;
186 }
187
188 if(do_data_recovery > 0){
189 ✗ z0 = (double*)calloc(size_z,sizeof(double));
190 ✗ z1 = (double*)calloc(size_z,sizeof(double));
191 ✗ assertStreamPrint(threadData,0!=z0,"calloc failed");
192 ✗ assertStreamPrint(threadData,0!=z1,"calloc failed");
193 }
194
195 ✗ functionODE_residual(data, threadData, x0, y0, z0);
196
197 ✗ x = data->localData[0]->realVars;
198
199 /* use actually value for xScaling */
200 ✗ for (i = 0; i < size_A; i++) {
201 ✗ modelica_real nominal = getNominalFromScalarIdx(data->simulationInfo, data->modelData, VAR_KIND_STATE, i);
202 ✗ xScaling[i] = fmax(nominal, fabs(x[i]));
203 }
204
205 /* solverData->f1 must be set outside this function based on x */
206 ✗ for(i = 0; i < size_A; i++) {
207 ✗ xsave = x[i];
208 ✗ delta_hh = delta_h * (fabs(xsave) + 1.0);
209 ✗ if (xsave + delta_hh >= getMaxFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, i)) {
210 ✗ delta_hh *= -1;
211 }
212 ✗ x[i] += delta_hh / xScaling[i];
213 /* Calculate scaled difference quotient */
214 ✗ delta_hh = 1. / delta_hh * xScaling[i];
215
216 ✗ functionODE_residual(data, threadData, x1, y1, z1);
217
218 ✗ for(j = 0; j < size_A; j++) {
219 ✗ k = i * size_A + j;
220 ✗ matrixA[k] = (x1[j] - x0[j]) * delta_hh;
221 }
222 ✗ for(j = 0; j < size_C; j++) {
223 ✗ k = i * size_C + j;
224 ✗ matrixC[k] = (y1[j] - y0[j]) * delta_hh;
225 }
226 ✗ if(do_data_recovery > 0){
227 ✗ for(j = 0; j < size_z; j++) {
228 ✗ k = i * size_z + j;
229 ✗ matrixCz[k] = (z1[j] - z0[j]) * delta_hh;
230 }
231 }
232 ✗ x[i] = xsave;
233 }
234
235 ✗ free(xScaling);
236 ✗ free(x0);
237 ✗ free(y0);
238 ✗ free(x1);
239 ✗ free(y1);
240 ✗ if(do_data_recovery > 0){
241 ✗ free(z0);
242 ✗ free(z1);
243 }
244
245 ✗ return 0;
246 }
247
248 ✗ int functionJacBD_num(DATA* data, threadData_t *threadData, double *matrixB, double *matrixD, double *matrixDz)
249 {
250 ✗ const double delta_h = numericalDifferentiationDeltaXlinearize;
251 double delta_hh;
252 double usave;
253 double* u;
254
255 int i,j,k;
256
257 int do_data_recovery = 0;
258 ✗ if(matrixDz){
259 do_data_recovery = 1;
260 }
261
262 ✗ int size_x = data->modelData->nStates;
263 ✗ int size_u = data->modelData->nInputVars;
264 ✗ int size_y = data->modelData->nOutputVars;
265 ✗ int size_z = data->modelData->nVariablesReal - 2*data->modelData->nStates;
266 ✗ double* x0 = (double*)calloc(size_x,sizeof(double));
267 ✗ double* y0 = (double*)calloc(size_y,sizeof(double));
268 ✗ double* x1 = (double*)calloc(size_x,sizeof(double));
269 ✗ double* y1 = (double*)calloc(size_y,sizeof(double));
270 double* z0 = 0;
271 double* z1 = 0;
272
273 ✗ assertStreamPrint(threadData,0!=x0,"calloc failed");
274 ✗ assertStreamPrint(threadData,0!=y0,"calloc failed");
275 ✗ assertStreamPrint(threadData,0!=x1,"calloc failed");
276 ✗ assertStreamPrint(threadData,0!=y1,"calloc failed");
277
278 ✗ if(do_data_recovery > 0){
279 ✗ z0 = (double*)calloc(size_z,sizeof(double));
280 ✗ z1 = (double*)calloc(size_z,sizeof(double));
281 ✗ assertStreamPrint(threadData,0!=z0,"calloc failed");
282 ✗ assertStreamPrint(threadData,0!=z1,"calloc failed");
283 }
284
285 ✗ functionODE_residual(data, threadData, x0, y0, z0);
286
287 ✗ u = data->simulationInfo->inputVars;
288
289 /* solverData->f1 must be set outside this function based on x */
290 ✗ for(i = 0; i < size_u; i++) {
291 ✗ usave = u[i];
292 ✗ delta_hh = delta_h * (fabs(usave) + 1.0);
293 ✗ u[i] += delta_hh;
294 ✗ delta_hh = 1. / delta_hh;
295
296 ✗ functionODE_residual(data, threadData, x1, y1, z1);
297
298 ✗ for(j = 0; j < size_x; j++) {
299 ✗ k = i * size_x + j;
300 ✗ matrixB[k] = (x1[j] - x0[j]) * delta_hh;
301 }
302 ✗ for(j = 0; j < size_y; j++) {
303 ✗ k = i * size_y + j;
304 ✗ matrixD[k] = (y1[j] - y0[j]) * delta_hh;
305 }
306 ✗ if(do_data_recovery > 0){
307 ✗ for(j = 0; j < size_z; j++) {
308 ✗ k = i * size_z + j;
309 ✗ matrixDz[k] = (z1[j] - z0[j]) * delta_hh;
310 }
311 }
312 ✗ u[i] = usave;
313 }
314
315 ✗ free(x0);
316 ✗ free(y0);
317 ✗ free(x1);
318 ✗ free(y1);
319 ✗ if(do_data_recovery > 0){
320 ✗ free(z0);
321 ✗ free(z1);
322 }
323
324 ✗ return 0;
325 }
326
327
328 /* Calculate the jacobian matrix by analytical finite difference */
329 ✗ int functionJacA(DATA* data, threadData_t *threadData, double* jac){
330
331 ✗ const int index = data->callback->INDEX_JAC_A;
332 ✗ JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
333 unsigned int i,j,k;
334 k = 0;
335 ✗ if (jacobian->constantEqns != NULL) {
336 ✗ jacobian->constantEqns(data, threadData, jacobian, NULL);
337 }
338
339 ✗ for(i=0; i < jacobian->sizeCols; i++)
340 {
341 ✗ jacobian->seedVars[i] = 1.0;
342 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_JAC))
343 {
344 printf("Caluculate one col:\n");
345 ✗ for(j=0; j < jacobian->sizeCols;j++)
346 {
347 ✗ infoStreamPrint(OMC_LOG_JAC,0,"seed: jacobian->seedVars[%d]= %f",j,jacobian->seedVars[j]);
348 }
349 }
350
351 ✗ data->callback->functionJacA_column(data, threadData, jacobian, NULL);
352
353 ✗ for(j = 0; j < jacobian->sizeRows; j++)
354 {
355 ✗ jac[k++] = jacobian->resultVars[j];
356 ✗ infoStreamPrint(OMC_LOG_JAC,0,"write in jac[%d]-[%d,%d]=%g from row[%d]=%g",k-1,i,j,jac[k-1],i,jacobian->resultVars[j]);
357 }
358
359 ✗ jacobian->seedVars[i] = 0.0;
360 }
361 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_JAC))
362 {
363 ✗ infoStreamPrint(OMC_LOG_JAC,0,"Print jac:");
364 ✗ for(i=0; i < jacobian->sizeRows;i++)
365 {
366 ✗ for(j=0; j < jacobian->sizeCols;j++) {
367 ✗ printf("% .5e ",jac[i+j*jacobian->sizeCols]);
368 }
369 printf("\n");
370 }
371 }
372
373 ✗ return 0;
374 }
375 ✗ int functionJacB(DATA* data, threadData_t *threadData, double* jac){
376
377 ✗ const int index = data->callback->INDEX_JAC_B;
378 ✗ JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
379
380 unsigned int i,j,k;
381 k = 0;
382 ✗ if (jacobian->constantEqns != NULL) {
383 ✗ jacobian->constantEqns(data, threadData, jacobian, NULL);
384 }
385
386 ✗ for(i=0; i < jacobian->sizeCols; i++)
387 {
388 ✗ jacobian->seedVars[i] = 1.0;
389 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_JAC))
390 {
391 printf("Caluculate one col:\n");
392 ✗ for(j=0; j < jacobian->sizeCols;j++)
393 {
394 ✗ infoStreamPrint(OMC_LOG_JAC,0,"seed: jacobian->seedVars[%d]= %f",j,jacobian->seedVars[j]);
395 }
396 }
397
398 ✗ data->callback->functionJacB_column(data, threadData, jacobian, NULL);
399
400 ✗ for(j = 0; j < jacobian->sizeRows; j++)
401 {
402 ✗ jac[k++] = jacobian->resultVars[j];
403 ✗ infoStreamPrint(OMC_LOG_JAC,0,"write in jac[%d]-[%d,%d]=%g from row[%d]=%g",k-1,i,j,jac[k-1],i,jacobian->resultVars[j]);
404 }
405
406 ✗ jacobian->seedVars[i] = 0.0;
407 }
408 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_JAC))
409 {
410 ✗ infoStreamPrint(OMC_LOG_JAC, 0, "Print jac:");
411 ✗ for(i=0; i < jacobian->sizeRows;i++)
412 {
413 ✗ for(j=0; j < jacobian->sizeCols;j++)
414 ✗ printf("% .5e ",jac[i+j*jacobian->sizeCols]);
415 printf("\n");
416 }
417 }
418
419 ✗ return 0;
420 }
421 ✗ int functionJacC(DATA* data, threadData_t *threadData, double* jac){
422
423 ✗ const int index = data->callback->INDEX_JAC_C;
424 ✗ JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
425 unsigned int i,j,k;
426 k = 0;
427 ✗ if (jacobian->constantEqns != NULL) {
428 ✗ jacobian->constantEqns(data, threadData, jacobian, NULL);
429 }
430
431 ✗ for(i=0; i < jacobian->sizeCols; i++)
432 {
433 ✗ jacobian->seedVars[i] = 1.0;
434 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_JAC))
435 {
436 printf("Caluculate one col:\n");
437 ✗ for(j=0; j < jacobian->sizeCols;j++)
438 ✗ infoStreamPrint(OMC_LOG_JAC,0,"seed: jacobian->seedVars[%d]= %f",j,jacobian->seedVars[j]);
439 }
440
441 ✗ data->callback->functionJacC_column(data, threadData, jacobian, NULL);
442
443 ✗ for(j = 0; j < jacobian->sizeRows; j++)
444 {
445 ✗ jac[k++] = jacobian->resultVars[j];
446 ✗ infoStreamPrint(OMC_LOG_JAC,0,"write in jac[%d]-[%d,%d]=%g from row[%d]=%g",k-1,i,j,jac[k-1],i,jacobian->resultVars[j]);
447 }
448
449 ✗ jacobian->seedVars[i] = 0.0;
450 }
451 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_JAC))
452 {
453 ✗ infoStreamPrint(OMC_LOG_JAC, 0, "Print jac:");
454 ✗ for(i=0; i < jacobian->sizeRows;i++)
455 {
456 ✗ for(j=0; j < jacobian->sizeCols;j++)
457 ✗ printf("% .5e ",jac[i+j*jacobian->sizeCols]);
458 printf("\n");
459 }
460 }
461
462 ✗ return 0;
463 }
464 ✗ int functionJacD(DATA* data, threadData_t *threadData, double* jac){
465
466 ✗ const int index = data->callback->INDEX_JAC_D;
467 ✗ JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]);
468 unsigned int i,j,k;
469 k = 0;
470 ✗ if (jacobian->constantEqns != NULL) {
471 ✗ jacobian->constantEqns(data, threadData, jacobian, NULL);
472 }
473
474 ✗ for(i=0; i < jacobian->sizeCols; i++)
475 {
476 ✗ jacobian->seedVars[i] = 1.0;
477 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_JAC))
478 {
479 printf("Caluculate one col:\n");
480 ✗ for(j=0; j < jacobian->sizeCols;j++) {
481 ✗ infoStreamPrint(OMC_LOG_JAC,0,"seed: jacobian->seedVars[%d]= %f",j,jacobian->seedVars[j]);
482 }
483 }
484
485 ✗ data->callback->functionJacD_column(data, threadData, jacobian, NULL);
486
487 ✗ for(j = 0; j < jacobian->sizeRows; j++)
488 {
489 ✗ jac[k++] = jacobian->resultVars[j];
490 ✗ infoStreamPrint(OMC_LOG_JAC,0, "write in jac[%d]-[%d,%d]=%g from row[%d]=%g",k-1,i,j,jac[k-1],i,jacobian->resultVars[j]);
491 }
492
493 ✗ jacobian->seedVars[i] = 0.0;
494 }
495 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_JAC))
496 {
497 ✗ infoStreamPrint(OMC_LOG_JAC, 0, "Print jac:");
498 ✗ for(i=0; i < jacobian->sizeRows;i++)
499 {
500 ✗ for(j=0; j < jacobian->sizeCols;j++)
501 ✗ printf("% .5e ",jac[i+j*jacobian->sizeCols]);
502 printf("\n");
503 }
504 }
505
506 ✗ return 0;
507 }
508
509
510
511 ✗ int linearize(DATA* data, threadData_t *threadData)
512 {
513 /* Check if data recovery is requested */
514 ✗ int do_data_recovery = omc_flag[FLAG_L_DATA_RECOVERY] ? 1 : 0;
515
516 /* init linearization sizes */
517 ✗ int size_A = data->modelData->nStates;
518 ✗ int size_Inputs = data->modelData->nInputVars;
519 ✗ int size_Outputs = data->modelData->nOutputVars;
520 ✗ int size_z = data->modelData->nVariablesReal - 2*data->modelData->nStates;
521 ✗ double* matrixA = (double*)calloc(size_A*size_A,sizeof(double));
522 ✗ double* matrixB = (double*)calloc(size_A*size_Inputs,sizeof(double));
523 ✗ double* matrixC = (double*)calloc(size_Outputs*size_A,sizeof(double));
524 ✗ double* matrixD = (double*)calloc(size_Outputs*size_Inputs,sizeof(double));
525 double* matrixCz = 0;
526 double* matrixDz = 0;
527 string strA, strB, strC, strD, strCz, strDz, strX, strU, strZ0, filename, ext;
528
529 ✗ assertStreamPrint(threadData,0!=matrixA,"calloc failed");
530 ✗ assertStreamPrint(threadData,0!=matrixB,"calloc failed");
531 ✗ assertStreamPrint(threadData,0!=matrixC,"calloc failed");
532 ✗ assertStreamPrint(threadData,0!=matrixD,"calloc failed");
533
534 ✗ if(do_data_recovery > 0){
535 ✗ matrixCz = (double*)calloc(size_z*size_A,sizeof(double));
536 ✗ matrixDz = (double*)calloc(size_z*size_Inputs,sizeof(double));
537 ✗ assertStreamPrint(threadData,0!=matrixCz,"calloc failed");
538 ✗ assertStreamPrint(threadData,0!=matrixDz,"calloc failed");
539 }
540
541 /* Need to do this before changing anything so that we get a proper z0 */
542 ✗ if(do_data_recovery > 0){
543 ✗ if(size_z){
544 ✗ strZ0 = "{" + array2string(&data->localData[0]->realVars[2*size_A], 1, size_z, data) + "}";
545 }else{
546 strZ0 = "zeros(0)";
547 }
548 }
549
550 /* Can currently only extract data recovery matrices Cz and Dz numerically, so we do this first if necessary */
551 ✗ if(do_data_recovery > 0 || data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sizeTmpVars == 0){
552 /* Calculate numeric Jacobian */
553 ✗ if(functionJacAC_num(data, threadData, matrixA, matrixC, matrixCz))
554 {
555 ✗ throwStreamPrint(threadData, "Error, can not get Matrix A or C ");
556 return 1;
557 }
558 ✗ if(functionJacBD_num(data, threadData, matrixB, matrixD, matrixDz))
559 {
560 ✗ throwStreamPrint(threadData, "Error, can not get Matrix B or D ");
561 return 1;
562 }
563 }
564
565 /* Check if symbolic Jacobian available, if it is then use it (overwriting A,B,C,D if also doing data recovery) */
566 ✗ if (data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sizeTmpVars > 0){
567 /* Retrieve symbolic Jacobian */
568 /* Determine Matrix A */
569 JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]);
570 ✗ if(!data->callback->initialAnalyticJacobianA(data, threadData, jacobian)){
571 ✗ assertStreamPrint(threadData,0==functionJacA(data, threadData, matrixA),"Error, can not get Matrix A ");
572 }
573
574 /* Determine Matrix B */
575 ✗ jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_B]);
576 ✗ if(!data->callback->initialAnalyticJacobianB(data, threadData, jacobian)){
577 ✗ assertStreamPrint(threadData,0==functionJacB(data, threadData, matrixB),"Error, can not get Matrix B ");
578 }
579
580 /* Determine Matrix C */
581 ✗ jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_C]);
582 ✗ if(!data->callback->initialAnalyticJacobianC(data, threadData, jacobian)){
583 ✗ assertStreamPrint(threadData,0==functionJacC(data, threadData, matrixC),"Error, can not get Matrix C ");
584 }
585
586 /* Determine Matrix D */
587 ✗ jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_D]);
588 ✗ if(!data->callback->initialAnalyticJacobianD(data, threadData, jacobian)){
589 ✗ assertStreamPrint(threadData,0==functionJacD(data, threadData, matrixD),"Error, can not get Matrix D ");
590 }
591 }
592 ✗ if (data->modelData->linearizationDumpLanguage != OMC_LINEARIZE_DUMP_LANGUAGE_PYTHON)
593 {
594
595 ✗ strA = array2string(matrixA, size_A, size_A, data);
596 ✗ strB = array2string(matrixB, size_A, size_Inputs, data);
597 ✗ strC = array2string(matrixC, size_Outputs, size_A, data);
598 ✗ strD = array2string(matrixD, size_Outputs, size_Inputs, data);
599 ✗ if (do_data_recovery > 0)
600 {
601 ✗ strCz = array2string(matrixCz, size_z, size_A, data);
602 ✗ strDz = array2string(matrixDz, size_z, size_Inputs, data);
603 }
604
605 // The empty array {} is not valid modelica, so we need to put something
606 // inside the curly braces for x0 and u0. {for i in in 1:0} will create an
607 // empty array if needed.
608 ✗ if (size_A)
609 {
610 // fix dummping julia vector braces
611 ✗ if (data->modelData->linearizationDumpLanguage == OMC_LINEARIZE_DUMP_LANGUAGE_JULIA)
612 ✗ strX = "[" + array2string(data->localData[0]->realVars, 1, size_A, data) + "]";
613 else
614 ✗ strX = "{" + array2string(data->localData[0]->realVars, 1, size_A, data) + "}";
615 }
616 else
617 {
618 strX = "zeros(0)";
619 }
620
621 ✗ if (size_Inputs)
622 {
623 // fix dummping julia vector braces
624 ✗ if (data->modelData->linearizationDumpLanguage == OMC_LINEARIZE_DUMP_LANGUAGE_JULIA)
625 ✗ strU = "[" + array2string(data->simulationInfo->inputVars, 1, size_Inputs, data) + "]";
626 else
627 ✗ strU = "{" + array2string(data->simulationInfo->inputVars, 1, size_Inputs, data) + "}";
628 }
629 else
630 {
631 strU = "zeros(0)";
632 }
633 }
634 else
635 {
636 // convert the matrices to Python format
637 //infoStreamPrint(OMC_LOG_STDOUT, 0, "Python selected");
638 ✗ strA = array2PythonString(matrixA, size_A, size_A);
639 ✗ strB = array2PythonString(matrixB, size_A, size_Inputs);
640 ✗ strC = array2PythonString(matrixC, size_Outputs, size_A);
641 ✗ strD = array2PythonString(matrixD, size_Outputs, size_Inputs);
642 ✗ if (do_data_recovery > 0)
643 {
644 ✗ strCz = array2PythonString(matrixCz, size_z, size_A);
645 ✗ strDz = array2PythonString(matrixDz, size_z, size_Inputs);
646 }
647 // strA = "[[-2.887152375617477, -1.62655852935388], [-2.380918056675567, -2.388394731625707]]";
648 //infoStreamPrint(OMC_LOG_STDOUT, 0, strA.c_str());
649 ✗ if (size_A)
650 ✗ strX = "[" + array2string(data->localData[0]->realVars, 1, size_A, data) + "]";
651 else
652 strX = "[0]";
653
654 ✗ if (size_Inputs)
655 ✗ strU = "[" + array2string(data->simulationInfo->inputVars, 1, size_Inputs, data) + "]";
656 else
657 strU = "[0]";
658 }
659
660 ✗ free(matrixA);
661 ✗ free(matrixB);
662 ✗ free(matrixC);
663 ✗ free(matrixD);
664 ✗ if(do_data_recovery > 0){
665 ✗ free(matrixCz);
666 ✗ free(matrixDz);
667 }
668 ✗ switch(data->modelData->linearizationDumpLanguage){
669 case OMC_LINEARIZE_DUMP_LANGUAGE_MODELICA: ext = ".mo"; break;
670 case OMC_LINEARIZE_DUMP_LANGUAGE_MATLAB: ext = ".m"; break;
671 case OMC_LINEARIZE_DUMP_LANGUAGE_JULIA: ext = ".jl"; break;
672 case OMC_LINEARIZE_DUMP_LANGUAGE_PYTHON: ext = ".py"; break;
673 }
674 /* ticket #5927: Don't use the model name to prevent bad names for certain languages. */
675 ✗ if (omc_flag[FLAG_OUTPUT_PATH]) {
676 ✗ filename = string(omc_flagValue[FLAG_OUTPUT_PATH]) + "/" + "linearized_model" + string(ext);
677 } else {
678 ✗ filename = "linearized_model" + string(ext);
679 }
680
681 ✗ FILE *fout = omc_fopen(filename.c_str(),"wb");
682 ✗ assertStreamPrint(threadData,0!=fout,"Cannot open File %s",filename.c_str());
683
684 const char* frame = NULL;
685 ✗ if(do_data_recovery > 0){
686 ✗ frame = data->callback->linear_model_datarecovery_frame();
687 fprintf(fout, frame, strX.c_str(), strU.c_str(), strZ0.c_str(), strA.c_str(), strB.c_str(), strC.c_str(), strD.c_str(), strCz.c_str(), strDz.c_str());
688 }else{
689 ✗ frame = data->callback->linear_model_frame();
690 ✗ fprintf(fout, frame, strX.c_str(), strU.c_str(), strA.c_str(), strB.c_str(), strC.c_str(), strD.c_str(), (double) data->simulationInfo->stopTime);
691 }
692
693 ✗ fflush(fout);
694 ✗ fclose(fout);
695
696 ✗ if (0 == strcmp(frame, "")) {
697 ✗ errorStreamPrint(OMC_LOG_STDOUT, 0, "Linear model could not be created.");
698 } else {
699 ✗ if (data->modelData->runTestsuite) {
700 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Linear model is created.");
701 }
702 else {
703 ✗ if (omc_flag[FLAG_OUTPUT_PATH]) {
704 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Linear model is created at %s", filename.c_str());
705 } else {
706 char* cwd = getcwd(NULL, 0); /* call with NULL and 0 to allocate the buffer dynamically (no pathmax needed) */
707 ✗ if(!cwd) {
708 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Linear model %s is created, but getting the full path failed.", filename.c_str());
709 }
710 else {
711 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Linear model is created at %s/%s", cwd, filename.c_str());
712 ✗ free(cwd);
713 }
714 }
715 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "The output format can be changed with the command line option --linearizationDumpLanguage.");
716 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "The options are: --linearizationDumpLanguage=none, modelica, matlab, julia, python.");
717 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "In OMEdit Simulation Setup->Linearize->Target language for linearized model.");
718 }
719 }
720 return 0;
721 }
722
723 }
724