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 / 218
Functions: 0.0% 0 / 0 / 20
Branches: 0.0% 0 / 0 / 310
OMCompiler/SimulationRuntime/c/moo/strategies.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 "../simulation/options.h"
29 #include "../simulation/arrayIndex.h"
30
31 #include "evaluations.h"
32
33 #include "strategies.h"
34
35 namespace OpenModelica {
36
37 // TODO: make more initializer methods: bionic?!
38
39 // ==================== Helpers for Emit and Simulation ====================
40
41 // use this field when needing some data object inside OpenModelica
42 // callbacks / interfaces but no nice void* field exists
43 static void *_global_reference_data_field = nullptr;
44
45 // set the global pointer
46 ✗ void set_global_reference_data(void *reference_data) {
47 ✗ assert(!_global_reference_data_field);
48 ✗ _global_reference_data_field = reference_data;
49 ✗ }
50
51 // get the global pointer
52 ✗ void* get_global_reference_data() {
53 ✗ assert(_global_reference_data_field);
54 ✗ return _global_reference_data_field;
55 }
56
57 // clear the global pointer
58 ✗ void clear_global_reference_data() {
59 ✗ _global_reference_data_field = nullptr;
60 ✗ }
61
62 ✗ static int control_trajectory_input_function(DATA* data, threadData_t* threadData) {
63 ✗ AuxiliaryControls* aux_controls = static_cast<AuxiliaryControls*>(get_global_reference_data());
64
65 ✗ const ControlTrajectory& controls = aux_controls->controls;
66 ✗ InfoGDOP& info = aux_controls->info;
67 f64* u_interpolation = aux_controls->u_interpolation.raw();
68
69 // important use data here not info.data
70 // for some reason a new data object is created for each solve
71 // but not important for now
72 ✗ f64 time = data->localData[0]->timeValue;
73
74 ✗ controls.interpolate_at(time, u_interpolation);
75 ✗ set_inputs(info, u_interpolation);
76
77 ✗ return 0;
78 }
79
80 ✗ static void trajectory_xut_emit(simulation_result* sim_result, DATA* data, threadData_t *threadData)
81 {
82 ✗ AuxiliaryTrajectory* aux = static_cast<AuxiliaryTrajectory*>(sim_result->storage); // exploit void* field
83 ✗ InfoGDOP& info = aux->info;
84 ✗ Trajectory& trajectory = aux->trajectory;
85 ✗ SOLVER_INFO* solver_info = aux->solver_info;
86
87 ✗ for (int x_idx = 0; x_idx < info.x_size; x_idx++) {
88 ✗ trajectory.x[x_idx].push_back(data->localData[0]->realVars[x_idx]);
89 }
90
91 ✗ for (int u_idx = 0; u_idx < info.u_size; u_idx++) {
92 ✗ int u = info.u_indices_real_vars[u_idx];
93 ✗ trajectory.u[u_idx].push_back(data->localData[0]->realVars[u]);
94 }
95
96 ✗ trajectory.t.push_back(solver_info->currentTime);
97 ✗ }
98
99 ✗ static void trajectory_p_emit(simulation_result* sim_result, DATA* data, threadData_t *threadData)
100 {
101 // TODO: PARAMETERS
102 ✗ AuxiliaryTrajectory* aux = static_cast<AuxiliaryTrajectory*>(sim_result->storage);
103 ✗ InfoGDOP& info = aux->info;
104 ✗ Trajectory& trajectory = aux->trajectory;
105
106 ✗ for (int p_idx = 0; p_idx < info.p_size; p_idx++) {
107 ✗ trajectory.p.push_back(data->localData[0]->realVars[0 /* parameter index */]);
108 }
109 ✗ }
110
111 // sets state and control initial values
112 // from e.g. initial equations / parameters
113 ✗ void initialize_model(InfoGDOP& info) {
114 ✗ externalInputallocate(info.data);
115 ✗ initializeModel(info.data, info.threadData, "", "", info.t0);
116 ✗ }
117
118 // at least free for externalInputallocate();
119 ✗ void free_model(InfoGDOP& info) {
120 ✗ externalInputFree(info.data);
121 ✗ }
122
123 // ==================== Emit to MAT file ====================
124
125 ✗ MatEmitter::MatEmitter(InfoGDOP& info) : info(info) {}
126
127 ✗ int MatEmitter::operator()(const PrimalDualTrajectory& trajectory) {
128 ✗ DATA* data = info.data;
129 ✗ threadData_t* threadData = info.threadData;
130
131 // TODO: this is placed poorly here -> maybe move the entry point from generated code deeper into the runtime
132 ✗ const char *result_file = omc_flagValue[FLAG_R];
133 std::string result_file_cstr;
134 ✗ if (result_file) {
135 ✗ data->modelData->resultFileName = GC_strdup(result_file);
136 ✗ } else if (omc_flag[FLAG_OUTPUT_PATH]) { /* read the output path from the command line (if any) */
137 ✗ if (0 > GC_asprintf(&result_file, "%s/%s_res.%s", omc_flagValue[FLAG_OUTPUT_PATH], data->modelData->modelFilePrefix, data->simulationInfo->outputFormat)) {
138 ✗ throwStreamPrint(NULL, "simulation_runtime.c: Error: can not allocate memory.");
139 }
140 ✗ data->modelData->resultFileName = GC_strdup(result_file);
141 } else {
142 ✗ result_file_cstr = std::string(data->modelData->modelFilePrefix) + std::string("_res.") + data->simulationInfo->outputFormat;
143 ✗ data->modelData->resultFileName = GC_strdup(result_file_cstr.c_str());
144 }
145
146 const auto& primals = trajectory.primals;
147
148 ✗ data->simulationInfo->numSteps = primals->t.size();
149 ✗ initializeResultData(data, threadData, 0);
150 ✗ sim_result.writeParameterData(&sim_result, data, threadData);
151
152 // allocate contiguous array for xu
153 ✗ FixedVector<f64> xu(primals->x.size() + primals->u.size());
154 ✗ for (size_t i = 0; i < primals->t.size(); i++) {
155 // move trajectory data in contiguous array
156 ✗ for (size_t x_index = 0; x_index < primals->x.size(); x_index++) {
157 ✗ xu[x_index] = primals->x[x_index][i];
158 }
159 ✗ for (size_t u_index = 0; u_index < primals->u.size(); u_index++) {
160 ✗ xu[primals->x.size() + u_index] = primals->u[u_index][i];
161 }
162
163 // evaluate all algebraic variables
164 ✗ set_time(info, primals->t[i]);
165 ✗ set_states_inputs(info, xu.raw());
166 ✗ eval_current_point_dae(info);
167
168 // emit point
169 ✗ sim_result.emit(&sim_result, data, threadData);
170 }
171
172 ✗ deinitializeResultData(data, threadData);
173
174 ✗ return 0;
175 }
176
177 // ==================== Constant Initialization ====================
178
179 ✗ ConstantInitialization::ConstantInitialization(InfoGDOP& info)
180 ✗ : info(info) {}
181
182 ✗ std::unique_ptr<PrimalDualTrajectory> ConstantInitialization::operator()(const GDOP::GDOP& gdop) {
183 ✗ DATA* data = info.data;
184
185 ✗ std::vector<f64> t = {info.t0, info.tf};
186 std::vector<std::vector<f64>> x_guess;
187 std::vector<std::vector<f64>> u_guess;
188 std::vector<f64> p;
189 ✗ InterpolationMethod interpolation = InterpolationMethod::LINEAR;
190
191 ✗ for (int x = 0; x < info.x_size; x++) {
192 ✗ if (data->modelData->realVarsData[x].dimension.numberOfDimensions > 0) {
193 ✗ Log::error("Support for array variables not yet implemented!");
194 ✗ abort();
195 }
196 ✗ modelica_real* start = (modelica_real *)data->modelData->realVarsData[x].attribute.start.data;
197 ✗ x_guess.push_back({start[0], start[0]});
198 }
199
200 ✗ for (int u : info.u_indices_real_vars) {
201 ✗ if (data->modelData->realVarsData[u].dimension.numberOfDimensions > 0) {
202 ✗ Log::error("Support for array variables not yet implemented!");
203 ✗ abort();
204 }
205 ✗ modelica_real* start = (modelica_real *)data->modelData->realVarsData[u].attribute.start.data;
206 ✗ u_guess.push_back({start[0], start[0]});
207 }
208
209 // TODO: PARAMETERS add p
210
211 ✗ return std::make_unique<PrimalDualTrajectory>(std::make_unique<Trajectory>(t, x_guess, u_guess, p, interpolation));
212 ✗ }
213
214 // ==================== Simulation ====================
215
216 ✗ Simulation::Simulation(InfoGDOP& info, SOLVER_METHOD solver)
217 ✗ : info(info), solver(solver) {}
218
219 ✗ std::unique_ptr<Trajectory> Simulation::operator()(const ControlTrajectory& controls, const FixedVector<f64>& parameters,
220 int num_steps, f64 start_time, f64 stop_time, f64* x_start_values) {
221 ✗ DATA* data = info.data;
222 ✗ threadData_t* threadData = info.threadData;
223 SOLVER_INFO solver_info;
224 ✗ SIMULATION_INFO *simInfo = data->simulationInfo;
225
226 ✗ solver_info.solverMethod = solver;
227 ✗ simInfo->numSteps = num_steps;
228 ✗ simInfo->startTime = start_time;
229 ✗ simInfo->stopTime = stop_time;
230 ✗ simInfo->stepSize = (stop_time - start_time) / static_cast<f64>(num_steps);
231 ✗ simInfo->useStopTime = 1;
232
233 // allocate and reserve trajectory vectors
234 std::vector<f64> t;
235 ✗ t.reserve(num_steps + 1);
236
237 ✗ std::vector<std::vector<f64>> x_sim(info.x_size);
238 ✗ for (auto& v : x_sim) v.reserve(num_steps + 1);
239
240 ✗ std::vector<std::vector<f64>> u_sim(info.u_size);
241 ✗ for (auto& v : u_sim) v.reserve(num_steps + 1);
242
243 ✗ std::vector<f64> p_sim(info.p_size);
244
245 // create Trajectory object
246 ✗ auto trajectory = std::make_unique<Trajectory>(Trajectory{t, x_sim, u_sim, p_sim, InterpolationMethod::LINEAR, nullptr});
247
248 // auxiliary data (passed as void* in storage member of sim_result)
249 ✗ auto aux_trajectory = std::make_unique<AuxiliaryTrajectory>(AuxiliaryTrajectory{*trajectory, info, &solver_info});
250
251 // define global sim_result
252 ✗ sim_result.filename = nullptr;
253 ✗ sim_result.numpoints = 0;
254 ✗ sim_result.cpuTime = 0;
255 ✗ sim_result.storage = aux_trajectory.get();
256 ✗ sim_result.emit = trajectory_xut_emit;
257 ✗ sim_result.init = nullptr;
258 ✗ sim_result.writeParameterData = trajectory_p_emit;
259 ✗ sim_result.free = nullptr;
260
261 // init simulation
262 ✗ initializeSolverData(data, threadData, &solver_info);
263 ✗ setZCtol(fmin(simInfo->stepSize, simInfo->tolerance));
264 ✗ initialize_model(info); // TODO: is this needed? we pass x0 after all, maybe call this when getting x0 from the model?!
265 ✗ data->real_time_sync.enabled = FALSE;
266
267 // create an auxiliary object (stored in global void*)
268 // since the input_function interface offers no additional argument
269 ✗ FixedVector<f64> u_interpolation_buffer = FixedVector<f64>(info.u_size);
270 ✗ auto aux_controls = AuxiliaryControls{controls, info, u_interpolation_buffer};
271 ✗ set_global_reference_data(&aux_controls);
272
273 // set the new input function
274 ✗ auto generated_input_function = data->callback->input_function;
275 ✗ data->callback->input_function = control_trajectory_input_function;
276
277 // set states and controls for start time
278 // emit for time = start_time
279 ✗ controls.interpolate_at(start_time, u_interpolation_buffer.raw());
280 ✗ set_inputs(info, u_interpolation_buffer.raw());
281 ✗ set_states(info, x_start_values);
282 ✗ eval_current_point_dae(info);
283 ✗ trajectory_xut_emit(&sim_result, data, threadData);
284
285 // ensure realVars stay consistent across ring buffer rotation (prefixedName_performSimulation line 491ff):
286 // copy current slot (localData[0]) into the upcoming slots (localData[1], localData[2])
287 // after rotateRingBuffer() + lookupRingBuffer(), the new "current slot"
288 // (localData[0]) will already contain our enforced values.
289 // necessary for simulation steps, as otherwise start values (t = 0) would be present in these slots!
290 ✗ memcpy(data->localData[1]->realVars, data->localData[0]->realVars, sizeof(modelica_real) * data->modelData->nVariablesReal);
291 ✗ memcpy(data->localData[2]->realVars, data->localData[1]->realVars, sizeof(modelica_real) * data->modelData->nVariablesReal);
292
293 // simulation with custom emit
294 ✗ data->callback->performSimulation(data, threadData, &solver_info);
295
296 // reset to previous input function (from generated code)
297 ✗ data->callback->input_function = generated_input_function;
298
299 // set global aux data to nullptr
300 ✗ clear_global_reference_data();
301
302 // free allocated memory
303 ✗ free_model(info); // free for initialize_model() (at least partial)
304
305 // Attention: this deletes the A Jacobian also!!
306 ✗ freeSolverData(data, &solver_info); // free for initializeSolverData
307
308 ✗ return trajectory;
309 ✗ }
310
311 // ==================== Simulation Step ====================
312
313 ✗ SimulationStep::SimulationStep(std::shared_ptr<Simulation> simulation) : simulation(simulation) {}
314
315 ✗ void SimulationStep::activate(const ControlTrajectory& controls_, const FixedVector<f64>& parameters_) {
316 ✗ controls = &controls_;
317 ✗ parameters = &parameters_;
318 ✗ }
319
320 ✗ void SimulationStep::reset() {
321 ✗ controls = nullptr;
322 ✗ parameters = nullptr;
323 ✗ }
324
325 ✗ std::unique_ptr<Trajectory> SimulationStep::operator()(f64* x_start_values, f64 start_time, f64 stop_time) {
326 ✗ if (!controls || !parameters) {
327 ✗ Log::error("OpenModelica::SimulationStep has not been activated.");
328 ✗ abort();
329 }
330 const int num_steps = 1;
331 ✗ return (*simulation)(*controls, *parameters, num_steps, start_time, stop_time, x_start_values);
332 }
333
334 // ==================== Nominal Scaling Factory ====================
335
336 ✗ std::shared_ptr<NLP::Scaling> NominalScalingFactory::operator()(const GDOP::GDOP& gdop) {
337 // x, g, f of the NLP { min f(x) s.t. g_l <= g(x) <= g_l }
338 ✗ auto x_nominal = FixedVector<f64>(gdop.get_number_vars());
339 ✗ auto g_nominal = FixedVector<f64>(gdop.get_number_constraints());
340 ✗ f64 f_nominal = 1;
341
342 // get problem sizes
343 ✗ auto x_size = info.x_size;
344 ✗ auto u_size = info.u_size;
345 ✗ auto xu_size = info.xu_size;
346 ✗ auto f_size = info.f_size;
347 ✗ auto g_size = info.g_size;
348 ✗ auto r_size = info.r_size;
349 ✗ auto fg_size = f_size + g_size;
350
351 ✗ auto has_mayer = gdop.get_problem().pc->has_mayer;
352 ✗ auto has_lagrange = gdop.get_problem().pc->has_lagrange;
353
354 ✗ if (has_mayer && has_lagrange) {
355 ✗ const modelica_real nominal_mayer = getNominalFromScalarIdx(info.data->simulationInfo, info.data->modelData, VAR_KIND_VARIABLE, info.index_mayer_real_vars);
356 ✗ const modelica_real nominal_lagrange = getNominalFromScalarIdx(info.data->simulationInfo, info.data->modelData, VAR_KIND_VARIABLE, info.index_lagrange_real_vars);
357
358 ✗ f_nominal = (nominal_mayer + nominal_lagrange) / 2;
359 }
360 ✗ else if (has_lagrange) {
361 ✗ f_nominal = getNominalFromScalarIdx(info.data->simulationInfo, info.data->modelData, VAR_KIND_VARIABLE, info.index_lagrange_real_vars);
362 }
363 ✗ else if (has_mayer) {
364 ✗ f_nominal = getNominalFromScalarIdx(info.data->simulationInfo, info.data->modelData, VAR_KIND_VARIABLE, info.index_mayer_real_vars);
365 }
366
367 // (x, u)_(t_node)
368 ✗ for (int node = 0; node < 1 + gdop.get_mesh().node_count; node++) {
369 ✗ for (int x = 0; x < x_size; x++) {
370 ✗ x_nominal[node * xu_size + x] = getNominalFromScalarIdx(info.data->simulationInfo, info.data->modelData, VAR_KIND_VARIABLE, x);
371 }
372
373 ✗ for (int u = 0; u < u_size; u++) {
374 ✗ int u_real_vars = info.u_indices_real_vars[u];
375 ✗ x_nominal[node * xu_size + x_size + u] = getNominalFromScalarIdx(info.data->simulationInfo, info.data->modelData, VAR_KIND_VARIABLE, u_real_vars);
376 }
377 }
378
379 ✗ for (int node = 0; node < gdop.get_mesh().node_count; node++) {
380 ✗ for (int f = 0; f < f_size; f++) {
381 ✗ g_nominal[node * fg_size + f] = x_nominal[f]; // reuse x nominal for dynamic for now!
382 }
383
384 ✗ for (int g = 0; g < g_size; g++) {
385 ✗ g_nominal[f_size + node * fg_size + g] = getNominalFromScalarIdx(info.data->simulationInfo, info.data->modelData, VAR_KIND_VARIABLE, info.index_g_real_vars + g);
386 }
387 }
388
389 ✗ for (int r = 0; r < r_size; r++) {
390 ✗ g_nominal[gdop.get_off_fg_total() + r] = getNominalFromScalarIdx(info.data->simulationInfo, info.data->modelData, VAR_KIND_VARIABLE, info.index_r_real_vars + r);
391 }
392
393 // artificial constraints are O(u)
394 ✗ for (int u = 0; u < info.u_size; u++) {
395 ✗ int u_real_vars = info.u_indices_real_vars[u];
396 ✗ g_nominal[gdop.get_off_fgr_total() + u] = getNominalFromScalarIdx(info.data->simulationInfo, info.data->modelData, VAR_KIND_VARIABLE, u_real_vars);
397 }
398
399 ✗ return std::make_shared<NLP::NominalScaling>(std::move(x_nominal), std::move(g_nominal), f_nominal);
400 }
401
402 // default strategies for OpenModelica
403 ✗ GDOP::Strategies default_strategies(InfoGDOP& info, GDOP::Problem& problem, bool use_moo_simulation) {
404 ✗ GDOP::Strategies strategies;
405
406 // TODO: do add simulation_tolerance factor here?
407 ✗ FixedVector<f64> verifier_tolerances(info.x_size);
408 ✗ for (int x = 0; x < info.x_size; x++) {
409 ✗ verifier_tolerances[x] = 1e-4 * getNominalFromScalarIdx(info.data->simulationInfo, info.data->modelData, VAR_KIND_VARIABLE, x);
410 }
411
412 auto scaling_factory = std::make_shared<NominalScalingFactory>(info);
413 ✗ auto emitter = std::make_shared<MatEmitter>(MatEmitter(info));
414 ✗ auto const_initialization_strategy = std::make_shared<ConstantInitialization>(ConstantInitialization(info));
415
416 ✗ std::shared_ptr<GDOP::Simulation> simulation_strategy;
417 ✗ std::shared_ptr<GDOP::SimulationStep> simulation_step_strategy;
418
419 ✗ if (!use_moo_simulation) {
420 ✗ auto tmp_simulation_strategy = std::make_shared<Simulation>(info, info.user_ode_solver);
421 ✗ simulation_step_strategy = std::make_shared<SimulationStep>(tmp_simulation_strategy);
422 simulation_strategy = tmp_simulation_strategy;
423 }
424 else {
425 ✗ simulation_strategy = std::make_shared<GDOP::RadauIntegratorSimulation>(*problem.dynamics);
426 ✗ simulation_step_strategy = std::make_shared<GDOP::RadauIntegratorSimulationStep>(*problem.dynamics);
427 }
428
429 ✗ auto simulation_initialization_strategy = std::make_shared<GDOP::SimulationInitialization>(GDOP::SimulationInitialization(const_initialization_strategy,
430 simulation_strategy));
431 ✗ auto verifier = std::make_shared<GDOP::SimulationVerifier>(GDOP::SimulationVerifier(simulation_strategy,
432 Linalg::Norm::NORM_INF,
433 std::move(verifier_tolerances)));
434
435 ✗ if (std::string(omc_flagValue[FLAG_IPOPT_INIT] ? omc_flagValue[FLAG_IPOPT_INIT] : "") == "CONST")
436 {
437 strategies.initialization = const_initialization_strategy;
438 }
439 else {
440 strategies.initialization = simulation_initialization_strategy;
441 }
442
443 strategies.simulation = simulation_strategy;
444 strategies.simulation_step = simulation_step_strategy;
445 ✗ strategies.mesh_refinement = std::make_shared<GDOP::L2BoundaryNorm>(info.l2bn_phase_one_iterations, info.l2bn_phase_two_iterations, info.l2bn_phase_two_level);
446 ✗ strategies.interpolation = std::make_shared<GDOP::PolynomialInterpolation>();
447 strategies.emitter = emitter;
448 strategies.verifier = verifier;
449 strategies.scaling_factory = scaling_factory;
450 ✗ strategies.refined_initialization = std::make_shared<GDOP::InterpolationRefinedInitialization>(
451 ✗ GDOP::InterpolationRefinedInitialization(strategies.interpolation, true, true, true));
452 ✗ return strategies;
453 ✗ };
454
455 } // namespace OpenModelica
456