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 / 176
Functions: 0.0% 0 / 0 / 5
Branches: 0.0% 0 / 0 / 96
OMCompiler/SimulationRuntime/c/optimization/DataManagement/InitialGuess.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 /*! InitialGuess.c
29 */
30
31 #include "../OptimizerData.h"
32 #include "../OptimizerLocalFunction.h"
33
34 #include "../../util/omc_file.h"
35
36 #include "simulation/arrayIndex.h"
37 #include "simulation/options.h"
38 #include "simulation/results/simulation_result.h"
39 #include "simulation/solver/dassl.h"
40 #include "simulation/solver/external_input.h"
41 #include "simulation/solver/initialization/initialization.h"
42 #include "simulation/solver/model_help.h"
43
44
45 static int initial_guess_ipopt_cflag(OptData *optData, char* cflags);
46 static inline void smallIntSolverStep(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo, const double tstop);
47 static short initial_guess_ipopt_sim(OptData *optData, SOLVER_INFO* solverInfo, const short o);
48 static inline void init_ipopt_data(OptData *optData, const short o);
49
50 /*!
51 * create initial guess
52 * author: Vitalij Ruge
53 **/
54 ✗ void initial_guess_optimizer(OptData *optData, SOLVER_INFO* solverInfo){
55
56 char *cflags;
57 int opt = 1;
58 int i, j;
59 char buffer[4096];
60 ✗ const int nu = optData->dim.nu;
61 ✗ optData->dim.iter = 0;
62 ✗ optData->ipop.csvOstep = (char*)(omc_flagValue[FLAG_CSV_OSTEP]);
63 ✗ optData->ipop.debugeJ = (char*)omc_flagValue[FLAG_OPTDEBUGEJAC];
64
65 ✗ optData->pFile = omc_fopen("optimizeInput.csv", "wt");
66
67 fprintf(optData->pFile, "%s ", "time");
68 ✗ for(i=0; i < nu; ++i){
69 ✗ sprintf(buffer, "%s", optData->dim.inputName[i]);
70 ✗ fprintf(optData->pFile, "%s ", buffer);
71 }
72 ✗ fprintf(optData->pFile, "%s", "\n");
73
74 ✗ cflags = (char*)omc_flagValue[FLAG_IPOPT_INIT];
75
76 ✗ if(cflags){
77 ✗ opt = initial_guess_ipopt_cflag(optData, cflags);
78 }
79
80 ✗ if(opt > 0)
81 ✗ opt = initial_guess_ipopt_sim(optData, solverInfo, opt);
82
83 ✗ init_ipopt_data(optData, opt);
84 ✗ }
85
86
87 /*!
88 * create initial guess dasslColorSymJac
89 * author: Vitalij Ruge
90 **/
91 ✗ static short initial_guess_ipopt_sim(OptData *optData, SOLVER_INFO* solverInfo, const short o)
92 {
93 double *u0;
94 int i,j,k,l;
95 modelica_real ***v;
96 long double tol;
97 short printGuess, op=1;
98
99 ✗ const int nx = optData->dim.nx;
100 ✗ const int nu = optData->dim.nu;
101 ✗ const int np = optData->dim.np;
102 ✗ const int nsi = optData->dim.nsi;
103 ✗ const int nReal = optData->dim.nReal;
104 ✗ char *cflags = (char*)omc_flagValue[FLAG_IIF];
105
106 ✗ DATA* data = optData->data;
107 ✗ threadData_t *threadData = optData->threadData;
108 ✗ SIMULATION_INFO *sInfo = data->simulationInfo;
109 const char *solverMethod = NULL;
110
111 ✗ if(!data->simulationInfo->external_input.active){
112 ✗ externalInputallocate(data);
113 }
114
115 /* Initial DASSL solver */
116 ✗ DASSL_DATA* dasslData = (DASSL_DATA*) malloc(sizeof(DASSL_DATA));
117 ✗ tol = data->simulationInfo->tolerance;
118 ✗ data->simulationInfo->tolerance = fmin(fmax(tol,1e-8),1e-3);
119
120 ✗ infoStreamPrint(OMC_LOG_SOLVER, 0, "Initial Guess: Initializing DASSL");
121 /* Borrowed for the duration of the guess; sInfo owns the name it came with. */
122 ✗ solverMethod = sInfo->solverMethod;
123 ✗ sInfo->solverMethod = "dassl";
124 ✗ solverInfo->solverMethod = S_DASSL;
125 ✗ dassl_initial(data, threadData, solverInfo, dasslData);
126 ✗ solverInfo->solverMethod = S_OPTIMIZATION;
127 ✗ solverInfo->solverData = dasslData;
128
129 ✗ u0 = optData->bounds.u0;
130 ✗ v = optData->v;
131
132 ✗ if(!data->simulationInfo->external_input.active)
133 ✗ for(i = 0; i< nu;++i)
134 ✗ data->simulationInfo->inputVars[i] = u0[i]/*optData->bounds.scalF[i + nx]*/;
135
136 ✗ printGuess = (short)(OMC_ACTIVE_STREAM(OMC_LOG_INIT) && !OMC_ACTIVE_STREAM(OMC_LOG_SOLVER));
137
138 ✗ if((double)data->simulationInfo->startTime < optData->time.t0){
139 double t = data->simulationInfo->startTime;
140
141 ✗ FILE * pFile = optData->pFile;
142 fprintf(pFile, "%lf ",(double)t);
143 ✗ for(i = 0; i < nu; ++i){
144 ✗ fprintf(pFile, "%lf ", (float)data->simulationInfo->inputVars[i]);
145 }
146 fprintf(pFile, "%s", "\n");
147 if(1){
148 printf("\nPreSim");
149 printf("\n========================================================\n");
150 ✗ printf("\ndone: time[%i] = %g",0,(double)data->simulationInfo->startTime);
151 }
152 ✗ while(t < optData->time.t0){
153 ✗ externalInputUpdate(data);
154 ✗ smallIntSolverStep(data, threadData, solverInfo, fmin(t += optData->time.dt[0], optData->time.t0));
155 ✗ printf("\ndone: time[%i] = %g",0,(double)data->localData[0]->timeValue);
156 ✗ sim_result.emit(&sim_result,data,threadData);
157 ✗ fprintf(pFile, "%lf ",(double)data->localData[0]->timeValue);
158 ✗ for(i = 0; i < nu; ++i){
159 ✗ fprintf(pFile, "%lf ", (float)data->simulationInfo->inputVars[i]);
160 }
161 fprintf(pFile, "%s", "\n");
162 }
163 ✗ copy_initial_values(optData, data);
164
165 if(1){
166 printf("\n--------------------------------------------------------");
167 printf("\nfinished: PreSim");
168 printf("\n========================================================\n");
169 }
170 }
171
172 ✗ if(o == 2 && cflags && strcmp(cflags, ""))
173 op = 2;
174
175 ✗ if(printGuess ){
176 printf("\nInitial Guess");
177 printf("\n========================================================\n");
178 ✗ printf("\ndone: time[%i] = %g",0,(double)optData->time.t0);
179 }
180
181 ✗ for(i = 0, k=1; i < nsi; ++i){
182 ✗ for(j = 0; j < np; ++j, ++k){
183 ✗ externalInputUpdate(data);
184 ✗ if(op==1)
185 ✗ smallIntSolverStep(data, threadData, solverInfo, (double)optData->time.t[i][j]);
186 else{
187 ✗ rotateRingBuffer(data->simulationData, 1);
188 ✗ lookupRingBuffer(data->simulationData, (void**) data->localData);
189 ✗ importStartValues(data, threadData, cflags, (double)optData->time.t[i][j]);
190 ✗ for(l=0; l<nReal; ++l){
191 ✗ data->localData[0]->realVars[l] = getStartFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, l);
192 }
193 }
194
195 ✗ if(printGuess)
196 ✗ printf("\ndone: time[%i] = %g", k, (double)optData->time.t[i][j]);
197
198 ✗ memcpy(v[i][j], data->localData[0]->realVars, nReal*sizeof(double));
199 ✗ for(l = 0; l < nx; ++l){
200
201 ✗ if(((double) v[i][j][l] < (double)optData->bounds.vmin[l]*optData->bounds.vnom[l])
202 ✗ || (double) (v[i][j][l] > (double) optData->bounds.vmax[l]*optData->bounds.vnom[l])){
203 printf("\n********************************************\n");
204 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Initial guess failure at time %g",(double)optData->time.t[i][j]);
205 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "%g<= (%s=%g) <=%g",
206 ✗ (double)optData->bounds.vmin[l]*optData->bounds.vnom[l],
207 ✗ data->modelData->realVarsData[l].info.name,
208 ✗ (double)v[i][j][l],
209 ✗ (double)optData->bounds.vmax[l]*optData->bounds.vnom[l]);
210 printf("\n********************************************");
211 }
212 }
213 }
214 }
215
216 ✗ if(printGuess){
217 printf("\n--------------------------------------------------------");
218 printf("\nfinished: Initial Guess");
219 printf("\n========================================================\n");
220 }
221
222 ✗ dassl_deinitial(data, solverInfo->solverData);
223 ✗ solverInfo->solverData = (void*)optData;
224 ✗ sInfo->solverMethod = solverMethod;
225 ✗ data->simulationInfo->tolerance = tol;
226
227 ✗ externalInputFree(data);
228 ✗ return op;
229 }
230
231
232 /*!
233 * helper for initial_guess_optimizer (pick up clfag option)
234 * author: Vitalij Ruge
235 **/
236 ✗ static int initial_guess_ipopt_cflag(OptData *optData, char* cflags)
237 {
238 ✗ if(!strcmp(cflags,"const") || !strcmp(cflags,"CONST"))
239 {
240 int i, j;
241 ✗ const int nsi = optData->dim.nsi;
242 ✗ const int np = optData->dim.np;
243 ✗ const int nu = optData->dim.nu;
244 ✗ const int nReal = optData->dim.nReal;
245
246 ✗ for(i = 0; i< nu; ++i )
247 ✗ optData->data->simulationInfo->inputVars[i] = optData->bounds.u0[i];
248 ✗ for(i = 0; i < nsi; ++i){
249 ✗ for(j = 0; j < np; ++j){
250 ✗ memcpy(optData->v[i][j], optData->v0, nReal*sizeof(modelica_real));
251 }
252 }
253
254 ✗ infoStreamPrint(OMC_LOG_IPOPT, 0, "Using const trajectory as initial guess.");
255 ✗ return 0;
256 ✗ }else if(!strcmp(cflags,"sim") || !strcmp(cflags,"SIM")){
257
258 ✗ infoStreamPrint(OMC_LOG_IPOPT, 0, "Using simulation as initial guess.");
259 ✗ return 1;
260 ✗ }else if(!strcmp(cflags,"file") || !strcmp(cflags,"FILE")){
261 ✗ infoStreamPrint(OMC_LOG_STDOUT, 0, "Using values from file as initial guess.");
262 ✗ return 2;
263 }
264
265 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "not support ipopt_init=%s", cflags);
266 ✗ return 1;
267
268 }
269
270 /*!
271 * init ipopt data struct
272 * author: Vitalij Ruge
273 **/
274 ✗ static inline void init_ipopt_data(OptData *optData, const short op){
275 OptDataIpopt* ipop = &optData->ipop;
276 ✗ DATA * data = optData->data;
277 ✗ const int NV = optData->dim.NV;
278 ✗ const int NRes = optData->dim.NRes;
279 ✗ const int nsi = optData->dim.nsi;
280 ✗ const int np = optData->dim.np;
281 ✗ const int nv = optData->dim.nv;
282 ✗ const int nc = optData->dim.nc;
283 ✗ const int ncf = optData->dim.ncf;
284 ✗ const int nJ = optData->dim.nJ;
285 ✗ const int nx = optData->dim.nx;
286 ✗ const int nReal = optData->dim.nReal;
287 ✗ const int index_con = optData->dim.index_con;
288 ✗ const int index_conf = optData->dim.index_conf;
289
290 int i,j,l,shift;
291
292 ✗ ipop->vopt = malloc(NV*sizeof(double));
293 ✗ ipop->mult_x_L = calloc(NV, sizeof(double));
294 ✗ ipop->mult_x_U = calloc(NV, sizeof(double));
295
296 ✗ ipop->gmin = calloc(NRes, sizeof(double));
297 ✗ ipop->gmax = calloc(NRes, sizeof(double));
298 ✗ ipop->mult_g = calloc(NRes, sizeof(double));
299
300 /* An OPT_LOOP_INPUT takes the value of the variable that replaced it when the
301 * guess comes from a file; the model names the pair, this applies it. */
302 ✗ int *uIdx = (int*) malloc((nv-nx)*sizeof(int));
303 ✗ int *uLoop = (int*) malloc((nv-nx)*sizeof(int));
304 ✗ data->callback->getInputVarIndicesInOptimization(data, uIdx, uLoop);
305
306 ✗ for(i = 0, shift = 0; i < nsi; ++i){
307 ✗ for(j = 0; j < np; ++j, shift+=nv){
308 ✗ memcpy(data->localData[0]->realVars, optData->v[i][j], nReal*sizeof(double));
309 ✗ if(op == 2){
310 ✗ for(l = 0; l < nv-nx; ++l){
311 ✗ if(uLoop[l] >= 0){
312 ✗ data->localData[0]->realVars[uIdx[l]] = data->localData[0]->realVars[uLoop[l]];
313 }
314 }
315 }
316 ✗ optData->data->callback->setInputData(optData->data);
317 ✗ for(l = 0; l<nx; ++l){
318 ✗ ipop->vopt[l + shift] = optData->v[i][j][l]*optData->bounds.scalF[l];
319 }
320 ✗ for(;l<nv;++l){
321 ✗ ipop->vopt[l + shift] = data->simulationInfo->inputVars[l-nx] * optData->bounds.scalF[l];
322 }
323 }
324 }
325
326
327 ✗ l = NRes-ncf;
328 ✗ for(j = 0; j< nc; ++j){
329 ✗ for(i = nx; i < l; i += nJ){
330 ✗ ipop->gmin[i+j] = getMinFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, j + index_con);
331 ✗ ipop->gmax[i+j] = getMaxFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, j + index_con);
332 }
333 }
334
335 /*terminal constraint(s)*/
336 ✗ for(j = 0; j < ncf; ++j, ++i){
337 ✗ ipop->gmin[l+j] = getMinFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, j + index_conf);
338 ✗ ipop->gmax[l+j] = getMaxFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, j + index_conf);
339 }
340
341 ✗ free(uIdx);
342 ✗ free(uLoop);
343 ✗ }
344
345 ✗ static inline void smallIntSolverStep(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo, const double tstop){
346 long double a;
347 int iter;
348 int err;
349
350 ✗ solverInfo->currentTime = data->localData[0]->timeValue;
351 ✗ while(solverInfo->currentTime < tstop){
352 a = 1.0;
353 iter = 0;
354
355 ✗ rotateRingBuffer(data->simulationData, 1);
356 ✗ lookupRingBuffer(data->simulationData, (void**) data->localData);
357 do{
358 ✗ if(data->modelData->nStates < 1){
359 ✗ solverInfo->currentTime = tstop;
360 ✗ data->localData[0]->timeValue = tstop;
361 ✗ break;
362 }
363 ✗ solverInfo->currentStepSize = a*(tstop - solverInfo->currentTime);
364 ✗ err = dassl_step(data, threadData, solverInfo);
365 ✗ a *= 0.5;
366 ✗ if(++iter > 10){
367 printf("\n");
368 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "Initial guess failure at time %.12g", solverInfo->currentTime);
369 ✗ assert(0);
370 }
371 ✗ }while(err < 0);
372
373 ✗ data->callback->updateContinuousSystem(data, threadData);
374
375 }
376 ✗ }
377