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 / 192
Functions: 0.0% 0 / 1 / 8
Branches: 0.0% 0 / 0 / 172
OMCompiler/SimulationRuntime/c/optimization/eval_all/EvalL.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 /*! EvalL.c
29 */
30
31 #include "../OptimizerData.h"
32 #include "../OptimizerLocalFunction.h"
33 #include "../../simulation/solver/model_help.h"
34
35 static inline void calculate_hessian_matrix_numerical(double * v, const double * const lambda, const double objFactor , OptData *optData, const int i, const int j);
36 static inline void calculate_weighted_sum_with_lagrange_multiplicator_from_tensor(const int i, const int j, double * res, const modelica_boolean upC, OptData *optData);
37 static inline void calculate_hessian_matrix_numerical_last_time_intervall(double * v, const double * const lambda, const double objFactor, OptData *optData, const int i, const int j);
38 static inline void calculate_weighted_sum_with_lagrange_multiplicator_from_tensor_last_time_intervall(const int i, const int j, double * res, const modelica_boolean upC, const modelica_boolean upC2, OptData *optData);
39
40
41 static inline long double guess_step_size_for_numerical_differentiation(const long double v)
42 {
43 ✗ return 1e-5*fabsl(v) + 1e-8;
44 }
45
46
47 ✗ static void init_hessian_structure(int *iRow, int *iCol, OptData *optData){
48 int i, j, k, p, l, r, c;
49 ✗ const int nsi = optData->dim.nsi;
50 ✗ const int np = optData->dim.np;
51 ✗ const int np1 = np + 1;
52 ✗ const int nv = optData->dim.nv;
53
54 ✗ for(i = 0, r = 0, c = 0, k = 0; i + 1 < nsi; ++i){
55 ✗ for(p = 1; p < np1; ++p, r += nv, c += nv){
56 ✗ for(j = 0; j < nv; ++j){
57 ✗ for(l = 0; l < j+1; ++l){
58 ✗ if(optData->s.H0[j][l]){
59 ✗ iRow[k] = r + j;
60 ✗ iCol[k++] = c + l;
61 }
62 }
63 }
64 }
65 }
66
67 /* init_hessian_structure_last_time_intervall */
68 ✗ for(p = 1; p < np1; ++p, r += nv, c += nv){
69 ✗ for(j = 0; j< nv; ++j){
70 ✗ for(l = 0; l< j+1; ++l){
71 ✗ if(optData->s.H1[j][l] && np == p){
72 ✗ iRow[k] = r + j;
73 ✗ iCol[k++] = c + l;
74 ✗ }else if(optData->s.H0[j][l]){
75 ✗ iRow[k] = r + j;
76 ✗ iCol[k++] = c + l;
77 }
78 }
79 }
80 }
81 ✗ }
82
83 ✗ static void fill_hessian_values(double *vopt, double *lambda, double obj_factor, OptData *optData, double *values){
84 ✗ const int nJ = optData->dim.nJ;
85 const modelica_boolean do_update_cost_function = obj_factor != 0;
86 ✗ const modelica_boolean do_update_mayer_cost_function = do_update_cost_function && optData->s.mayer;
87 ✗ const modelica_boolean do_update_lagrange_cost_function = do_update_cost_function && optData->s.lagrange;
88 ✗ const int np = optData->dim.np;
89 ✗ const int np1 = np + 1;
90 ✗ const int nv = optData->dim.nv;
91 ✗ const int nsi = optData->dim.nsi;
92
93 int ii, p, i, j, k;
94 double * v;
95 double * la;
96
97 ✗ DATA * data = optData->data;
98
99 ✗ ++optData->dim.iter;
100 ✗ optData->dim.iter_updateHessian = 0;
101 ✗ copy_initial_values(optData, data);
102
103 ✗ for(ii = 0, k = 0, v = vopt, la = lambda; ii + 1 < nsi; ++ii){
104 ✗ for(p = 1; p < np1; ++p, v += nv, la += nJ){
105 ✗ calculate_hessian_matrix_numerical(v, la, obj_factor, optData, ii, p-1);
106 /*******************/
107 ✗ for(i = 0; i < nv; ++i){
108 ✗ for(j = 0; j < i + 1; ++j){
109 ✗ if(optData->s.H0[i][j]){
110 ✗ calculate_weighted_sum_with_lagrange_multiplicator_from_tensor(i, j, values + (k++), do_update_lagrange_cost_function, optData);
111 }
112 }
113 }
114 /*******************/
115 }
116 }
117
118 /* fill_hessian_values last time intervall */
119 ✗ for(p = 1; p < np1; ++p, v += nv, la += nJ){
120 ✗ calculate_hessian_matrix_numerical_last_time_intervall(v, la, obj_factor, optData, ii, p-1);
121 ✗ for(i = 0; i < nv; ++i){
122 ✗ for(j = 0; j < i + 1; ++j){
123 ✗ if(optData->s.H1[i][j] && np == p){
124 ✗ calculate_weighted_sum_with_lagrange_multiplicator_from_tensor_last_time_intervall(i, j, values + (k++), do_update_lagrange_cost_function, do_update_mayer_cost_function, optData);
125 ✗ }else if(optData->s.H0[i][j]){
126 ✗ calculate_weighted_sum_with_lagrange_multiplicator_from_tensor(i, j, values + (k++),do_update_lagrange_cost_function, optData);
127 }
128 }
129 }
130 }
131
132 ✗ }
133
134 /* eval hessian
135 */
136 ✗ Bool ipopt_h(int n, double *vopt, Bool new_x, double obj_factor, int m, double *lambda, Bool new_lambda,
137 int nele_hess, int *iRow, int *iCol, double *values, void* useData){
138
139 OptData *optData = (OptData*)useData;
140 ✗ modelica_boolean keepH = optData->dim.updateHessian > optData->dim.iter_updateHessian++;
141
142 ✗ if(values == NULL){
143 ✗ init_hessian_structure(iRow, iCol, optData);
144 ✗ }else if(keepH){
145 ✗ memcpy(values,optData->oldH,nele_hess*sizeof(double));
146 }else{
147 ✗ if(optData->ipop.csvOstep)
148 ✗ debugeSteps(optData, vopt, lambda);
149 ✗ fill_hessian_values(vopt,lambda, obj_factor, optData, values);
150 ✗ if(optData->dim.updateHessian > 0)
151 ✗ memcpy(optData->oldH, values, nele_hess*sizeof(double));
152 }
153
154
155 ✗ return TRUE;
156 }
157
158 /* numerical approximation
159 * hessian
160 */
161 ✗ static inline void calculate_hessian_matrix_numerical(double * v, const double * const lambda,
162 const double objFactor , OptData *optData, const int i, const int j){
163
164 ✗ const modelica_boolean la = optData->s.lagrange;
165 ✗ const modelica_boolean upCost = la && objFactor != 0;
166 ✗ DATA * data = optData->data;
167 ✗ threadData_t *threadData = optData->threadData;
168
169 ✗ const int nv = optData->dim.nv;
170 ✗ const int nx = optData->dim.nx;
171 ✗ const int nJ = optData->dim.nJ;
172 ✗ const modelica_real * const vmax = optData->bounds.vmax;
173 ✗ const modelica_real * const vnom = optData->bounds.vnom;
174
175 int ii,jj, l;
176 long double v_save, h;
177 modelica_real * realV[3];
178
179
180 ✗ for(l = 1; l<3; ++l){
181 ✗ realV[l] = data->localData[l]->realVars;
182 ✗ data->localData[l]->realVars = optData->v[i][j];
183 ✗ data->localData[l]->timeValue = (modelica_real) optData->time.t[i][j];
184 }
185 ✗ data->localData[0]->timeValue = (modelica_real) optData->time.t[i][j];
186
187 ✗ for(ii = 0; ii < nv; ++ii){
188 /********************/
189 ✗ v_save = (long double) v[ii];
190 h = (long double)guess_step_size_for_numerical_differentiation(v_save);
191 ✗ v[ii] += h;
192 ✗ if( v[ii] >= vmax[ii]){
193 ✗ h *= -1.0;
194 ✗ v[ii] = v_save + h;
195 }
196 /********************/
197 ✗ for(l = 0; l < nx; ++l)
198 ✗ data->localData[0]->realVars[l] = v[l]*vnom[l];
199 ✗ for(; l <nv; ++l)
200 ✗ data->simulationInfo->inputVars[l-nx] = (modelica_real) v[l]*vnom[l];
201 ✗ data->callback->input_function(data, threadData);
202 /*data->callback->functionDAE(data);*/
203 ✗ updateDiscreteSystem(data, threadData);
204 /********************/
205 ✗ diffSynColoredOptimizerSystem(optData, optData->tmpJ, i,j,2);
206 /********************/
207 ✗ v[ii] = (double)v_save;
208 /********************/
209 ✗ for(jj = 0; jj <ii+1; ++jj){
210 ✗ if(optData->s.H0[ii][jj]){
211 ✗ for(l = 0; l < nJ; ++l){
212 ✗ if(optData->s.Hg[l][ii][jj] && lambda[l] != 0)
213 ✗ optData->H[l][ii][jj] = (long double)(optData->tmpJ[l][jj] - optData->J[i][j][l][jj])*lambda[l]/h;
214 }
215 }
216 }
217 /********************/
218 ✗ if(upCost){
219 ✗ h = objFactor/h;
220 ✗ for(jj = 0; jj <ii+1; ++jj){
221 ✗ if(optData->s.Hl[ii][jj]){
222 ✗ optData->Hl[ii][jj] = (long double)(optData->tmpJ[nJ][jj] - optData->J[i][j][nJ][jj])*h;
223 }else{
224 ✗ optData->Hl[ii][jj] = 0.0;
225 }
226 }
227 }
228 /********************/
229 }
230
231 ✗ for(l = 1; l<3; ++l){
232 ✗ data->localData[l]->realVars = realV[l];
233 }
234
235 ✗ }
236
237
238 /* numerical approximation
239 * hessian
240 */
241 ✗ static inline void calculate_hessian_matrix_numerical_last_time_intervall(double * v, const double * const lambda,
242 const double objFactor, OptData *optData, const int i, const int j){
243
244 ✗ const modelica_boolean la = optData->s.lagrange;
245 ✗ const modelica_boolean ma = optData->s.mayer;
246 ✗ const modelica_boolean upCost = la && objFactor != 0;
247
248 ✗ const int nv = optData->dim.nv;
249 ✗ const int nx = optData->dim.nx;
250 ✗ const int np = optData->dim.np;
251 ✗ const int nJ = optData->dim.nJ;
252 ✗ const int nJ1 = optData->dim.nJ + 1;
253 ✗ const int nsi = optData->dim.nsi;
254 ✗ const int ncf = optData->dim.ncf;
255 ✗ const modelica_real * const vmax = optData->bounds.vmax;
256 ✗ const modelica_real * const vnom = optData->bounds.vnom;
257 ✗ const modelica_boolean upFinalCon = np == j + 1 && nsi == i +1;
258 ✗ const modelica_boolean upCost2 = upFinalCon && ma && objFactor != 0;
259 const short indexJ = (upCost2) ? 3 : 2;
260 int ii,jj, l,k;
261 long double v_save, h;
262 ✗ DATA * data = optData->data;
263 ✗ threadData_t *threadData = optData->threadData;
264
265 modelica_real * realV[3];
266
267 ✗ data->localData[0]->timeValue = (modelica_real) optData->time.t[i][j];
268 ✗ for(l = 1; l<3; ++l){
269 ✗ realV[l] = data->localData[l]->realVars;
270 ✗ data->localData[l]->realVars = optData->v[i][j];
271 ✗ data->localData[l]->timeValue = (modelica_real) optData->time.t[i][j];
272 }
273
274 ✗ for(ii = 0; ii < nv; ++ii){
275 /********************/
276 ✗ v_save = (long double) v[ii];
277 h = (long double)guess_step_size_for_numerical_differentiation(v_save);
278 ✗ v[ii] += h;
279 ✗ if(v[ii] > vmax[ii]){
280 ✗ h *= -1.0;
281 ✗ v[ii] = v_save + h;
282 }
283 /********************/
284 ✗ for(l = 0; l < nx; ++l)
285 ✗ data->localData[0]->realVars[l] = v[l]*vnom[l];
286 ✗ for(; l <nv; ++l)
287 ✗ data->simulationInfo->inputVars[l-nx] = (modelica_real) v[l]*vnom[l];
288
289 ✗ data->callback->input_function(data, threadData);
290 /*data->callback->functionDAE(data);*/
291 ✗ updateDiscreteSystem(data, threadData);
292 /********************/
293 ✗ diffSynColoredOptimizerSystem(optData, optData->tmpJ, i,j,indexJ);
294 /********************/
295 ✗ v[ii] = (double)v_save;
296 /********************/
297 ✗ for(jj = 0; jj <ii+1; ++jj){
298 ✗ if(optData->s.H0[ii][jj]){
299 ✗ for(l = 0; l < nJ; ++l){
300 ✗ if(optData->s.Hg[l][ii][jj])
301 ✗ optData->H[l][ii][jj] = (long double)(optData->tmpJ[l][jj] - optData->J[i][j][l][jj])*lambda[l]/h;
302 }
303 }
304 }
305 /********************/
306 ✗ if(upCost){
307 long double hh;
308 ✗ hh = objFactor/h;
309 ✗ for(jj = 0; jj <ii+1; ++jj){
310 ✗ if(optData->s.Hl[ii][jj]){
311 ✗ optData->Hl[ii][jj] = (long double)(optData->tmpJ[nJ][jj] - optData->J[i][j][nJ][jj])*hh;
312 }
313 }
314 }
315 /********************/
316 ✗ if(upCost2){
317 long double hh;
318 ✗ hh = objFactor/h;
319 ✗ for(jj = 0; jj <ii+1; ++jj){
320 ✗ if(optData->s.Hm[ii][jj]){
321 ✗ optData->Hm[ii][jj] = (long double)(optData->tmpJ[nJ1][jj] - optData->J[i][j][nJ1][jj])*hh;
322 }
323 }
324 }
325 /********************/
326 ✗ if(upFinalCon && ncf > 0){
327 ✗ diffSynColoredOptimizerSystemF(optData, optData->tmpJf);
328 ✗ for(jj = 0; jj <ii+1; ++jj){
329 ✗ if(optData->s.H0[ii][jj]){
330 ✗ for(l = 0; l < ncf; ++l){
331 ✗ if(optData->s.Hcf[l][ii][jj]){
332 ✗ optData->Hcf[l][ii][jj] = (long double)(optData->tmpJf[l][jj] - optData->Jf[l][jj])*lambda[nJ+l]/h;
333 }
334 }
335 }
336 }
337 }
338 /********************/
339 }
340
341 ✗ for(l = 1; l<3; ++l){
342 ✗ data->localData[l]->realVars = realV[l];
343 }
344
345 ✗ }
346
347
348 /* eval hessian for lagrange
349 */
350 ✗ static inline void calculate_weighted_sum_with_lagrange_multiplicator_from_tensor(const int i, const int j, double * res,
351 const modelica_boolean upC, OptData *optData){
352 ✗ const int nJ = optData->dim.nJ;
353
354 long double sum = 0.0;
355 int l;
356
357 ✗ for(l = 0; l< nJ; ++l){
358 ✗ if(optData->s.Hg[l][i][j])
359 ✗ sum += optData->H[l][i][j];
360 }
361
362 ✗ if(upC && optData->s.Hl[i][j])
363 ✗ sum += optData->Hl[i][j];
364
365 ✗ *res = (double) sum;
366
367 ✗ }
368
369 /* eval hessian for lagrange
370 */
371 ✗ static inline void calculate_weighted_sum_with_lagrange_multiplicator_from_tensor_last_time_intervall(const int i, const int j, double * res,
372 const modelica_boolean upC, const modelica_boolean upC2, OptData *optData){
373 ✗ const int nJ = optData->dim.nJ;
374 ✗ const int ncf = optData->dim.ncf;
375
376 long double sum = 0.0;
377 int l;
378
379 ✗ if(optData->s.H0[i][j]){
380 ✗ for(l = 0; l< nJ; ++l){
381 ✗ if(optData->s.Hg[l][i][j])
382 ✗ sum += optData->H[l][i][j];
383 }
384
385 ✗ if(upC && optData->s.Hl[i][j])
386 ✗ sum += optData->Hl[i][j];
387 }
388 ✗ for(l = 0; l< ncf; ++l){
389 ✗ if(optData->s.Hcf[l][i][j])
390 ✗ sum += optData->Hcf[l][i][j];
391 }
392 ✗ if(upC2 && optData->s.Hm[i][j])
393 ✗ sum += optData->Hm[i][j];
394
395 ✗ *res = (double) sum;
396
397 ✗ }
398
399