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 / 215
Functions: 0.0% 0 / 0 / 4
Branches: 0.0% 0 / 0 / 128
OMCompiler/SimulationRuntime/c/optimization/optimizer_main.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 /*! optimizer_main.c
29 * move model data in optimizer structure
30 */
31
32 #include "OptimizerData.h"
33 #include "OptimizerLocalFunction.h"
34 #include "simulation_data.h"
35 #include "simulation/options.h"
36
37 static inline int optimizationWithIpopt(OptData*optData);
38 static inline void freeOptimizerData(OptData*optData);
39 static const char* firstArrayVariable(const MODEL_DATA *modelData);
40
41 ✗ int runOptimizer(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo){
42 OptData *optData, optData_;
43 const char *arrayVariable;
44
45 /* The optimizer maps variables to optimization variables by scalar index. */
46 ✗ arrayVariable = firstArrayVariable(data->modelData);
47 ✗ if (arrayVariable != NULL) {
48 ✗ throwStreamPrint(threadData, "Optimization does not support array variables, but %s is an array. "
49 "Use --simCodeScalarize=true.", arrayVariable);
50 }
51
52 ✗ solverInfo->solverData = &optData_;
53 ✗ data->simulationInfo->noThrowDivZero = 1;
54
55 ✗ pickUpModelData(data, threadData, solverInfo);
56 ✗ optData = (OptData*) solverInfo->solverData;
57
58 ✗ initial_guess_optimizer(optData, solverInfo);
59 ✗ allocate_der_struct(&optData->s, &optData->dim ,data, optData);
60
61 ✗ const int res = optimizationWithIpopt(optData);
62 ✗ res2file(optData, solverInfo, optData->ipop.vopt);
63 ✗ freeOptimizerData(optData);
64 ✗ if(res == 0 /*Solve_Succeeded*/ || res == 1 /*Solved_To_Acceptable_Level*/)
65 ✗ return 0;
66 return -1;
67 }
68
69 /**
70 * @brief Find first array variable or parameter.
71 *
72 * @param modelData Model data.
73 * @return const char* Name of first array variable or parameter,
74 * NULL if all variables are scalars.
75 */
76 ✗ static const char* firstArrayVariable(const MODEL_DATA *modelData)
77 {
78 long i;
79
80 #define FIND_ARRAY(VARS, N) \
81 for (i = 0; i < (N); ++i) { \
82 if ((VARS)[i].dimension.numberOfDimensions > 0) { \
83 return (VARS)[i].info.name; \
84 } \
85 }
86
87 ✗ FIND_ARRAY(modelData->realVarsData, modelData->nVariablesRealArray)
88 ✗ FIND_ARRAY(modelData->integerVarsData, modelData->nVariablesIntegerArray)
89 ✗ FIND_ARRAY(modelData->booleanVarsData, modelData->nVariablesBooleanArray)
90 ✗ FIND_ARRAY(modelData->stringVarsData, modelData->nVariablesStringArray)
91 ✗ FIND_ARRAY(modelData->realParameterData, modelData->nParametersRealArray)
92 ✗ FIND_ARRAY(modelData->integerParameterData, modelData->nParametersIntegerArray)
93 ✗ FIND_ARRAY(modelData->booleanParameterData, modelData->nParametersBooleanArray)
94 ✗ FIND_ARRAY(modelData->stringParameterData, modelData->nParametersStringArray)
95
96 #undef FIND_ARRAY
97
98 return NULL;
99 }
100
101 /*!
102 * run optimization with ipopt
103 * author: Vitalij Ruge
104 **/
105 ✗ static inline int optimizationWithIpopt(OptData*optData){
106 IpoptProblem nlp = NULL;
107
108 ✗ const int NV = optData->dim.NV;
109 ✗ const int NRes = optData->dim.NRes;
110 ✗ const int nsi = optData->dim.nsi;
111 ✗ const int np = optData->dim.np;
112 ✗ const int nx = optData->dim.nx;
113 ✗ const int NJ = optData->dim.nJderx;
114 ✗ const int NJf = optData->dim.nJfderx;
115 ✗ const int nH0 = optData->dim.nH0_;
116 ✗ const int nH1 = optData->dim.nH1_;
117 ✗ const int njac = np*(NJ*nsi + nx*(np*nsi - 1)) + NJf;
118 ✗ const int nhess = (nsi*np-1)*nH0+nH1;
119
120 ✗ ipnumber * Vmin = optData->bounds.Vmin;
121 ✗ ipnumber * Vmax = optData->bounds.Vmax;
122 ✗ ipnumber * gmin = optData->ipop.gmin;
123 ✗ ipnumber * gmax = optData->ipop.gmax;
124 ✗ ipnumber * vopt = optData->ipop.vopt;
125 ✗ ipnumber * mult_g = optData->ipop.mult_g;
126 ✗ ipnumber * mult_x_L = optData->ipop.mult_x_L;
127 ✗ ipnumber * mult_x_U = optData->ipop.mult_x_U;
128 ipnumber obj;
129
130 char *cflags;
131 int max_iter = 5000;
132 int res = 0;
133
134 ✗ nlp = CreateIpoptProblem(NV, Vmin, Vmax,
135 NRes, gmin, gmax, njac, nhess, 0, &evalfF,
136 &evalfG, &evalfDiffF, &evalfDiffG, &ipopt_h);
137
138 /********************************************************************/
139 /******************* ipopt flags ************************/
140 /********************************************************************/
141
142 /*tol */
143 ✗ AddIpoptNumOption(nlp, "tol", optData->data->simulationInfo->tolerance);
144 ✗ AddIpoptStrOption(nlp, "evaluate_orig_obj_at_resto_trial", "yes");
145
146 /* print level */
147 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_IPOPT_FULL)){
148 ✗ AddIpoptIntOption(nlp, "print_level", 7);
149 ✗ }else if(OMC_ACTIVE_STREAM(OMC_LOG_IPOPT)){
150 ✗ AddIpoptIntOption(nlp, "print_level", 5);
151 ✗ }else if(OMC_ACTIVE_STREAM(OMC_LOG_STATS)){
152 ✗ AddIpoptIntOption(nlp, "print_level", 3);
153 }else {
154 ✗ AddIpoptIntOption(nlp, "print_level", 2);
155 }
156 ✗ AddIpoptIntOption(nlp, "file_print_level", 0);
157
158 /* derivative_test */
159 ✗ if(OMC_ACTIVE_STREAM(OMC_LOG_IPOPT_JAC) && OMC_ACTIVE_STREAM(OMC_LOG_IPOPT_HESSE)){
160 ✗ AddIpoptIntOption(nlp, "print_level", 4);
161 ✗ AddIpoptStrOption(nlp, "derivative_test", "second-order");
162 ✗ }else if(OMC_ACTIVE_STREAM(OMC_LOG_IPOPT_JAC)){
163 ✗ AddIpoptIntOption(nlp, "print_level", 4);
164 ✗ AddIpoptStrOption(nlp, "derivative_test", "first-order");
165 ✗ }else if(OMC_ACTIVE_STREAM(OMC_LOG_IPOPT_HESSE)){
166 ✗ AddIpoptIntOption(nlp, "print_level", 4);
167 ✗ AddIpoptStrOption(nlp, "derivative_test", "only-second-order");
168 }else{
169 ✗ AddIpoptStrOption(nlp, "derivative_test", "none");
170 }
171
172
173 ✗ cflags = (char*)omc_flagValue[FLAG_IPOPT_HESSE];
174 ✗ if(cflags){
175 ✗ if(!strcmp(cflags,"BFGS"))
176 ✗ AddIpoptStrOption(nlp, "hessian_approximation", "limited-memory");
177 ✗ else if(!strcmp(cflags,"const") || !strcmp(cflags,"CONST"))
178 ✗ AddIpoptStrOption(nlp, "hessian_constant", "yes");
179 ✗ else if(!(!strcmp(cflags,"num") || !strcmp(cflags,"NUM")))
180 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "not support ipopt_hesse=%s",cflags);
181 }
182
183 /*linear_solver e.g. mumps, MA27, MA57,...
184 * be sure HSL solver are installed if your try HSL solver*/
185 ✗ cflags = (char*)omc_flagValue[FLAG_LS_IPOPT];
186 ✗ if(cflags)
187 ✗ AddIpoptStrOption(nlp, "linear_solver", cflags);
188 ✗ AddIpoptNumOption(nlp,"mumps_pivtolmax",1e-5);
189
190
191 /* max iter */
192 ✗ cflags = (char*)omc_flagValue[FLAG_IPOPT_MAX_ITER];
193 ✗ if(cflags){
194 char buffer[100];
195 char c;
196 int index_e = -1, i = 0;
197 strcpy(buffer,cflags);
198
199 ✗ while(buffer[i] != '\0'){
200 ✗ if(buffer[i] == 'e'){
201 index_e = i;
202 break;
203 }
204 ✗ ++i;
205 }
206
207 ✗ if(index_e < 0){
208 max_iter = atoi(cflags);
209 ✗ if(max_iter >= 0)
210 ✗ AddIpoptIntOption(nlp, "max_iter", max_iter);
211 printf("\nmax_iter = %i",atoi(cflags));
212
213 }else{
214 ✗ max_iter = (atoi(cflags)*pow(10.0, (double)atoi(cflags+index_e+1)));
215 ✗ if(max_iter >= 0)
216 ✗ AddIpoptIntOption(nlp, "max_iter", (int)max_iter);
217 ✗ printf("\nmax_iter = (int) %i | (double) %g",(int)max_iter, atoi(cflags)*pow(10.0, (double)atoi(cflags+index_e+1)));
218 }
219 }else
220 ✗ AddIpoptIntOption(nlp, "max_iter", 5000);
221
222 /*heuristic optition */
223 {
224 int ws = 0;
225 ✗ cflags = (char*)omc_flagValue[FLAG_IPOPT_WARM_START];
226 ✗ if(cflags){
227 ws = atoi(cflags);
228 }
229
230 ✗ if(ws > 0){
231 ✗ double shift = pow(10,-1.0*ws);
232 ✗ AddIpoptNumOption(nlp,"mu_init",shift);
233 ✗ AddIpoptNumOption(nlp,"bound_mult_init_val",shift);
234 ✗ AddIpoptStrOption(nlp,"mu_strategy", "monotone");
235 ✗ AddIpoptNumOption(nlp,"bound_push", 1e-5);
236 ✗ AddIpoptNumOption(nlp,"bound_frac", 1e-5);
237 ✗ AddIpoptNumOption(nlp,"slack_bound_push", 1e-5);
238 ✗ AddIpoptNumOption(nlp,"constr_mult_init_max", 1e-5);
239 ✗ AddIpoptStrOption(nlp,"bound_mult_init_method","mu-based");
240 }else{
241 ✗ AddIpoptStrOption(nlp,"mu_strategy","adaptive");
242 ✗ AddIpoptStrOption(nlp,"bound_mult_init_method","constant");
243 }
244 ✗ AddIpoptStrOption(nlp,"fixed_variable_treatment","make_parameter");
245 ✗ AddIpoptStrOption(nlp,"dependency_detection_with_rhs","yes");
246 ✗ AddIpoptNumOption(nlp,"nu_init",1e-9);
247 ✗ AddIpoptNumOption(nlp,"eta_phi",1e-10);
248 }
249
250 /********************************************************************/
251
252
253 ✗ if(max_iter >=0){
254 ✗ optData->iter_ = 0.0;
255 ✗ optData->index = 1;
256 ✗ res = IpoptSolve(nlp, vopt, NULL, &obj, mult_g, mult_x_L, mult_x_U, (void*)optData);
257 }
258 ✗ if(res != 0 && !OMC_ACTIVE_STREAM(OMC_LOG_IPOPT))
259 ✗ warningStreamPrint(OMC_LOG_STDOUT, 0, "No optimal solution found!\nUse -lv=LOG_IPOPT for more information.");
260 ✗ FreeIpoptProblem(nlp);
261 ✗ return res;
262 }
263
264
265 ✗ static inline void freeOptimizerData(OptData*optData){
266 ✗ const int nsi = optData->dim.nsi;
267 ✗ const int np = optData->dim.np;
268 ✗ const int nv = optData->dim.nv;
269 ✗ const int nJ = optData->dim.nJ;
270
271 int i,j,k;
272
273 /*************************/
274 ✗ for(i=0; i < nsi; ++i)
275 ✗ free(optData->time.t[i]);
276 ✗ free(optData->time.t);
277 ✗ free(optData->time.dt);
278 /*************************/
279 ✗ free(optData->bounds.vmin);
280 ✗ free(optData->bounds.vmax);
281 ✗ free(optData->bounds.Vmin);
282 ✗ free(optData->bounds.Vmax);
283 ✗ free(optData->bounds.vnom);
284 ✗ free(optData->bounds.scalF);
285 ✗ for(i=0; i < nsi; ++i)
286 ✗ free(optData->bounds.scaldt[i]);
287 ✗ free(optData->bounds.scaldt);
288 ✗ for(i = 0; i < nsi; ++i){
289 ✗ free(optData->bounds.scalb[i]);
290 }
291 ✗ free(optData->bounds.scalb);
292 ✗ free(optData->bounds.u0);
293 /*************************/
294 ✗ free(optData->ipop.vopt);
295 ✗ free(optData->ipop.gmin);
296 ✗ free(optData->ipop.gmax);
297 ✗ free(optData->ipop.mult_g);
298 ✗ free(optData->ipop.mult_x_L);
299 ✗ free(optData->ipop.mult_x_U);
300 /*************************/
301
302 ✗ for(k = 0; k < nv; ++k){
303 ✗ free(optData->s.H0[k]);
304 ✗ free(optData->s.H1[k]);
305 ✗ free(optData->s.Hm[k]);
306 ✗ free(optData->s.Hl[k]);
307 }
308 ✗ free(optData->s.H0);
309 ✗ free(optData->s.H1);
310 ✗ free(optData->s.Hm);
311 ✗ free(optData->s.Hl);
312 ✗ for(j = 0; j < nJ; ++j){
313 ✗ for(k = 0; k < nv; ++k)
314 ✗ free(optData->s.Hg[j][k]);
315 ✗ free(optData->s.Hg[j]);
316 }
317 ✗ free(optData->s.Hg);
318 ✗ free(optData->s.lindex);
319 ✗ free(optData->s.seedVec);
320 ✗ free(optData->s.indexCon2);
321 ✗ free(optData->s.indexCon3);
322 ✗ free(optData->s.indexJ2);
323 ✗ free(optData->s.indexJ3);
324
325 /*************************/
326 ✗ for(i = 0; i < nsi; ++i){
327 ✗ for(j = 0; j < np; ++j){
328 ✗ free(optData->v[i][j]);
329 }
330 ✗ free(optData->v[i]);
331 }
332 ✗ free(optData->v);
333 ✗ free(optData->v0);
334 ✗ free(optData->sv0);
335 ✗ free(optData->b0);
336 ✗ free(optData->i0);
337 ✗ free(optData->b0Pre);
338 ✗ free(optData->i0Pre);
339 ✗ free(optData->v0Pre);
340 ✗ free(optData->rePre);
341 ✗ free(optData->re);
342 ✗ free(optData->storeR);
343
344 ✗ for(i = 0; i < nsi; ++i){
345 ✗ for(j = 0; j < np; ++j){
346 ✗ for(k = 0; k < nJ; ++k)
347 ✗ free(optData->J[i][j][k]);
348 ✗ free(optData->J[i][j]);
349 }
350 ✗ free(optData->J[i]);
351 }
352 ✗ free(optData->J);
353 ✗ for(k = 0; k < nJ; ++k)
354 ✗ free(optData->tmpJ[k]);
355 ✗ for(k = 0; k < nJ; ++k){
356 ✗ for(j = 0; j < nv; ++j){
357 ✗ free(optData->H[k][j]);
358 }
359 ✗ free(optData->H[k]);
360 }
361 ✗ free(optData->H);
362 ✗ for(j = 0; j < nv; ++j){
363 ✗ free(optData->Hl[j]);
364 ✗ free(optData->Hm[j]);
365 }
366 ✗ free(optData->Hl);
367 ✗ free(optData->Hm);
368 ✗ if(optData->dim.updateHessian > 0)
369 ✗ free(optData->oldH);
370
371 ✗ free(optData->dim.inputName);
372
373 ✗ for(k = 0; k < 2; ++k){
374 ✗ if(optData->s.matrix[2+k]){
375 ✗ for(i = 0; i< nsi; ++i){
376 ✗ for(j = 0; j< np; ++j){
377 ✗ free(optData->dim.analyticJacobians_tmpVars[k][i][j]);
378 }
379 ✗ free(optData->dim.analyticJacobians_tmpVars[k][i]);
380 }
381 ✗ free(optData->dim.analyticJacobians_tmpVars[k]);
382 }
383 }
384 ✗ free(optData->dim.analyticJacobians_tmpVars);
385
386 ✗ for(i = 0; i< optData->dim.nJ; ++i)
387 ✗ free(optData->s.JderCon[i]);
388 ✗ free(optData->s.JderCon);
389 ✗ free(optData->s.gradM);
390 ✗ free(optData->s.gradL);
391 ✗ }
392