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 / 181
Functions: 0.0% 0 / 0 / 14
Branches: 0.0% 0 / 0 / 150
OMCompiler/SimulationRuntime/c/moo/problem.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_data.h"
29 #include "simulation/arrayIndex.h"
30 #include "simulation/simulation_runtime.h"
31 #include "simulation/solver/external_input.h"
32 #include "simulation/solver/gbode_main.h"
33
34 #include <nlp/instances/gdop/gdop.h>
35
36 #include "evaluations.h"
37 #include "strategies.h"
38
39 #include "problem.h"
40
41 #define NUM_HES_FD_STEP 1e-8 // base step size for numerical Hessian perturbation
42 #define NUM_HES_DF_EXTR_STEPS 1 // number of extrapolation steps
43 #define NUM_HES_EXTR_DIV 2 // divisor for new step size in extrapolation
44
45
46 namespace OpenModelica {
47
48 ✗ FullSweep::FullSweep(GDOP::FullSweepLayout&& layout_lfg,
49 const GDOP::ProblemConstants& pc,
50 ✗ InfoGDOP& info)
51 ✗ : GDOP::FullSweep(std::move(layout_lfg), pc), info(info) {
52 ✗ }
53
54 ✗ void FullSweep::callback_eval(const f64* xu_nlp, const f64* p) {
55 ✗ set_parameters(info, p);
56 ✗ for (int i = 0; i < pc.mesh->intervals; i++) {
57 ✗ for (int j = 0; j < pc.mesh->nodes[i]; j++) {
58 ✗ const f64* xu_ij = get_xu_ij(xu_nlp, i, j);
59 ✗ f64* eval_buf_ij = get_eval_buffer(i, j);
60
61 ✗ set_states_inputs(info, xu_ij);
62 ✗ set_time(info, pc.mesh->t[i][j]);
63 ✗ eval_current_point_dae(info);
64 ✗ eval_lfg_write(info, eval_buf_ij);
65 }
66 }
67 ✗ }
68
69 ✗ void FullSweep::callback_jac(const f64* xu_nlp, const f64* p) {
70 ✗ set_parameters(info, p);
71 ✗ for (int i = 0; i < pc.mesh->intervals; i++) {
72 ✗ for (int j = 0; j < pc.mesh->nodes[i]; j++) {
73 ✗ const f64* xu_ij = get_xu_ij(xu_nlp, i, j);
74 ✗ f64* jac_buf_ij = get_jac_buffer(i, j);
75
76 ✗ set_states_inputs(info, xu_ij);
77 ✗ set_time(info, pc.mesh->t[i][j]);
78 ✗ eval_current_point_dae(info);
79 /* TODO: check if B matrix does hold additional ders */
80 ✗ jac_eval_write_as_csc(info, info.exc_jac->B.jacobian, jac_buf_ij);
81 }
82 }
83 ✗ }
84
85 ✗ void FullSweep::callback_hes(const f64* xu_nlp, const f64* p, const FixedField<f64, 2>& lagrange_factors, const f64* lambda) {
86 ✗ set_parameters(info, p);
87 ✗ for (int i = 0; i < pc.mesh->intervals; i++) {
88 ✗ for (int j = 0; j < pc.mesh->nodes[i]; j++) {
89 ✗ const f64* xu_ij = get_xu_ij(xu_nlp, i, j);
90 ✗ const f64* lambda_ij = get_lambda_ij(lambda, i, j);
91 ✗ f64* jac_buf_ij = get_jac_buffer(i, j);
92 ✗ f64* hes_buf_ij = get_hes_buffer(i, j);
93
94 ✗ set_states_inputs(info, xu_ij);
95 ✗ set_time(info, pc.mesh->t[i][j]);
96
97 /* TODO: check if B matrix does hold additional ders */
98
99 ✗ if (pc.has_lagrange) {
100 /* OpenModelica sorts the Functions as fLg, we have to
101 * use the workspace buffer from info.exc_hes and split the old lambda
102 * (lmbd_f_1, ..., lmbd_f_n, lmbd_g_1, ..., lmbd_g_m) in the middle as
103 * (lmbd_f_1, ..., lmbd_f_n, *L_factor*, lmbd_g_1, ..., lmbd_g_m) */
104 ✗ info.exc_hes->B.lambda[pc.x_size] = lagrange_factors[i][j];
105 ✗ for (int f = 0; f < pc.f_size; f++) {
106 ✗ info.exc_hes->B.lambda[f] = lambda_ij[f];
107 }
108 ✗ for (int g = 0; g < pc.g_size; g++) {
109 /* Lagrange offset */
110 ✗ info.exc_hes->B.lambda[pc.f_size + 1 + g] = lambda_ij[pc.f_size + g];
111 }
112 /* set wrapper lambda */
113 ✗ info.exc_hes->B.args.lambda = info.exc_hes->B.lambda.raw();
114 }
115 else {
116 /* set wrapper lambda, as Lagrange term isnt set and OM and MOO sortings are the same (f, g) */
117 ✗ info.exc_hes->B.args.lambda = lambda_ij;
118 }
119
120 /* set previous Jacobian *CSC* OpenModelica buffer */
121 ✗ info.exc_hes->B.args.jac_csc = jac_buf_ij;
122
123 /* call Hessian */
124 ✗ richardson_extrapolation(info.exc_hes->B.extr, hessian_fwd_differences_wrapper, &info.exc_hes->B.args,
125 NUM_HES_FD_STEP, NUM_HES_DF_EXTR_STEPS, NUM_HES_EXTR_DIV, 1, hes_buf_ij);
126 }
127 }
128 ✗ }
129
130 ✗ BoundarySweep::BoundarySweep(GDOP::BoundarySweepLayout&& layout_mr,
131 const GDOP::ProblemConstants& pc,
132 ✗ InfoGDOP& info)
133 ✗ : GDOP::BoundarySweep(std::move(layout_mr), pc), info(info) {}
134
135 ✗ void BoundarySweep::callback_eval(const f64* xu0_nlp, const f64* xuf_nlp, const f64* p, const f64 t0, const f64 tf) {
136 ✗ set_parameters(info, p);
137 ✗ set_states_inputs(info, xuf_nlp);
138 ✗ set_time(info, tf);
139 ✗ eval_current_point_dae(info);
140 ✗ eval_mr_write(info, get_eval_buffer());
141 ✗ }
142
143 ✗ void BoundarySweep::callback_jac(const f64* xu0_nlp, const f64* xuf_nlp, const f64* p, const f64 t0, const f64 tf) {
144 f64* jac_buf = get_jac_buffer();
145 ✗ set_parameters(info, p);
146 ✗ set_states_inputs(info, xuf_nlp);
147 ✗ set_time(info, tf);
148 ✗ eval_current_point_dae(info);
149 /* TODO: check if C matrix does hold additional ders */
150
151 /* derivative of mayer to jacbuffer[0] ... jac_buffer[exc_jac.D_coo.nnz_offset - 1] */
152 ✗ if (pc.has_mayer) {
153 ✗ jac_eval_write_first_row_as_csc(info, info.exc_jac->C.jacobian, info.exc_jac->C.buffer.raw(),
154 ✗ jac_buf, info.exc_jac->C.sparsity);
155 }
156
157 ✗ if (info.exc_jac->D.exists) {
158 /* TODO: check if D matrix does hold additional ders */
159 ✗ jac_eval_write_as_csc(info, info.exc_jac->D.jacobian, jac_buf + info.exc_jac->D.sparsity.nnz_offset);
160 }
161 ✗ }
162
163 ✗ void BoundarySweep::callback_hes(const f64* xu0_nlp, const f64* xuf_nlp, const f64* p, const f64 t0, const f64 tf, const f64 mayer_factor, const f64* lambda) {
164 ✗ set_parameters(info, p);
165 ✗ set_states_inputs(info, xuf_nlp);
166 ✗ set_time(info, tf);
167 fill_zero_hes_buffer();
168
169 f64* jac_buf = get_jac_buffer();
170 f64* hes_buffer = get_hes_buffer();
171
172 ✗ if (pc.has_mayer) {
173 // set all lambdas to 0, except mayer lambda
174 ✗ int index_mayer = pc.x_size + static_cast<int>(info.lagrange_exists);
175 ✗ info.exc_hes->C.lambda[index_mayer] = mayer_factor;
176 ✗ info.exc_hes->C.args.lambda = info.exc_hes->C.lambda.raw();
177
178 ✗ richardson_extrapolation(info.exc_hes->C.extr, hessian_fwd_differences_wrapper, &info.exc_hes->C.args,
179 NUM_HES_FD_STEP, NUM_HES_DF_EXTR_STEPS, NUM_HES_EXTR_DIV, 1, info.exc_hes->C.buffer.raw());
180 ✗ for (auto& [index_C, index_buffer] : info.exc_hes->C_to_Mr_buffer) {
181 ✗ hes_buffer[index_buffer] += info.exc_hes->C.buffer[index_C];
182 }
183 }
184
185 ✗ if (pc.r_size != 0) {
186 /* set duals and precomputed Jacobian D */
187 ✗ info.exc_hes->D.args.lambda = lambda;
188 ✗ info.exc_hes->D.args.jac_csc = jac_buf + info.exc_jac->D.sparsity.nnz_offset;
189
190 ✗ richardson_extrapolation(info.exc_hes->D.extr, hessian_fwd_differences_wrapper, &info.exc_hes->D.args,
191 NUM_HES_FD_STEP, NUM_HES_DF_EXTR_STEPS, NUM_HES_EXTR_DIV, 1, info.exc_hes->D.buffer.raw());
192 ✗ for (auto& [index_D, index_buffer] : info.exc_hes->D_to_Mr_buffer) {
193 ✗ hes_buffer[index_buffer] += info.exc_hes->D.buffer[index_D];
194 }
195 }
196 ✗ }
197
198 ✗ Dynamics::Dynamics(const GDOP::ProblemConstants& pc, InfoGDOP& info)
199 ✗ : GDOP::Dynamics(pc), info(info) {}
200
201 ✗ void Dynamics::allocate() {
202 // TODO: its unclear if other allocations are missing here. for now its working
203 ✗ JACOBIAN* jacobian = info.exc_jac->A.jacobian;
204 ✗ if (!jacobian->sparsePattern) {
205 ✗ info.data->callback->initialAnalyticJacobianA(info.data, info.threadData, jacobian);
206 ✗ allocated_ode_matrix = true;
207 }
208
209 ✗ SPARSE_PATTERN* sparsePattern = jacobian->sparsePattern;
210
211 ✗ jac_pattern = ::Simulation::Jacobian::sparse(
212 ::Simulation::JacobianFormat::CSC,
213 ✗ reinterpret_cast<int*>(sparsePattern->index) /* row indices */,
214 ✗ reinterpret_cast<int*>(sparsePattern->leadindex) /* col pointers */,
215 ✗ sparsePattern->nnz);
216 ✗ }
217
218 ✗ void Dynamics::free() {
219 ✗ if (allocated_ode_matrix) {
220 ✗ freeJacobian(info.exc_jac->A.jacobian);
221 ✗ allocated_ode_matrix = false;
222 }
223 ✗ }
224
225 ✗ void Dynamics::eval(const f64* x, const f64* u, const f64* p, f64 t, f64* f, void* user_data) {
226 ✗ set_parameters(info, p);
227 ✗ set_states(info, x);
228 ✗ set_inputs(info, u);
229 ✗ set_time(info, t);
230 ✗ eval_current_point_ode(info);
231 ✗ eval_ode_write(info, f);
232 ✗ }
233
234 ✗ void Dynamics::jac(const f64* x, const f64* u, const f64* p, f64 t, f64* dfdx, void* user_data) {
235 ✗ set_parameters(info, p);
236 ✗ set_states(info, x);
237 ✗ set_inputs(info, u);
238 ✗ set_time(info, t);
239 ✗ eval_write_ode_jacobian(info, dfdx);
240 ✗ }
241
242 ✗ GDOP::Problem create_gdop(InfoGDOP& info, Mesh& mesh) {
243 ✗ DATA* data = info.data;
244
245 // at first call init for all start values
246 ✗ initialize_model(info);
247
248 /* variable sizes */
249 ✗ info.x_size = data->modelData->nStates;
250 ✗ info.u_size = data->modelData->nInputVars;
251 ✗ info.xu_size = info.x_size + info.u_size;
252 ✗ info.p_size = 0; // TODO: PARAMETERS
253
254 /* variable bounds */
255 ✗ FixedVector<Bounds> x_bounds(info.x_size);
256 ✗ FixedVector<Bounds> u_bounds(info.u_size);
257 ✗ FixedVector<Bounds> p_bounds(info.p_size);
258
259 ✗ for (int x = 0; x < info.x_size; x++) {
260 ✗ x_bounds[x].lb = getMinFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_STATE, x);
261 ✗ x_bounds[x].ub = getMaxFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_STATE, x);
262 }
263
264 /* new generated function getInputVarIndices, just fills the index list of all optimizable inputs */
265 ✗ info.u_indices_real_vars = FixedVector<int>(info.u_size);
266 /* The loop-input table is the classic optimizer's; MOO does not use it. */
267 ✗ FixedVector<int> loopInputs(info.u_size);
268 ✗ data->callback->getInputVarIndicesInOptimization(data, info.u_indices_real_vars.raw(), loopInputs.raw());
269 ✗ for (int u = 0; u < info.u_size; u++) {
270 ✗ int u_index = info.u_indices_real_vars[u];
271 ✗ u_bounds[u].lb = getMinFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, u_index);
272 ✗ u_bounds[u].ub = getMaxFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, u_index);
273 }
274
275 /* constraint sizes */
276 ✗ info.f_size = info.x_size;
277 ✗ info.index_der_x_real_vars = info.x_size;
278 ✗ info.g_size = data->modelData->nOptimizeConstraints;
279 ✗ info.r_size = data->modelData->nOptimizeFinalConstraints; // TODO: add *generic boundary* constraints later also at t=t0
280
281 ✗ short der_index_mayer_realVars = -1;
282 ✗ short der_indices_lagrange_realVars[2] = {-1, -1};
283
284 /* this is really ugly IMO, fix this when ready for master! */
285 ✗ info.mayer_exists = (data->callback->mayer(data, &info.address_mayer_real_vars, &der_index_mayer_realVars) >= 0);
286 ✗ if (info.mayer_exists) {
287 ✗ info.index_mayer_real_vars = static_cast<int>(info.address_mayer_real_vars - data->localData[0]->realVars);
288 }
289
290 ✗ info.lagrange_exists = (data->callback->lagrange(data, &info.address_lagrange_real_vars, &der_indices_lagrange_realVars[0], &der_indices_lagrange_realVars[1]) >= 0);
291 ✗ if (info.lagrange_exists) {
292 ✗ info.index_lagrange_real_vars = static_cast<int>(info.address_lagrange_real_vars - data->localData[0]->realVars);
293 }
294
295 /* constraint bounds */
296 ✗ FixedVector<Bounds> g_bounds(info.g_size);
297 ✗ FixedVector<Bounds> r_bounds(info.r_size);
298
299 ✗ info.index_g_real_vars = data->modelData->nVariablesReal - (info.g_size + info.r_size);
300 ✗ info.index_r_real_vars = data->modelData->nVariablesReal - info.r_size;
301 ✗ for (int g = 0; g < info.g_size; g++) {
302 ✗ g_bounds[g].lb = getMinFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, info.index_g_real_vars + g);
303 ✗ g_bounds[g].ub = getMaxFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, info.index_g_real_vars + g);
304 }
305
306 ✗ for (int r = 0; r < info.r_size; r++) {
307 ✗ r_bounds[r].lb = getMinFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, info.index_r_real_vars + r);
308 ✗ r_bounds[r].ub = getMaxFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, info.index_r_real_vars + r);
309 }
310
311 /* for now we ignore xf fixed (need some steps in Backend to detect)
312 * and also ignore x0 non fixed, since too complicated
313 * => assume x(t_0) = x0 fixed, x(t_f) free to r constraint / maybe the old BE can do that already?!
314 * option: generate fixed final states individually */
315 ✗ FixedVector<std::optional<f64>> xu0_fixed(info.xu_size);
316 ✗ FixedVector<std::optional<f64>> xuf_fixed(info.xu_size);
317
318 /* fix time horizon for now to [t0, tf] */
319 ✗ std::array<Bounds, 2> T_bounds = { Bounds{ info.t0, info.t0 }, Bounds{ info.tf, info.tf } };
320 std::array<std::optional<f64>, 2> T_fixed = { info.t0, info.tf };
321
322 /* set *fixed* initial, final states */
323 ✗ for (int x = 0; x < info.x_size; x++) {
324 ✗ xu0_fixed[x] = getStartFromScalarIdx(data->simulationInfo, data->modelData, VAR_TYPE_REAL, VAR_KIND_VARIABLE, x);
325 }
326
327 /* u0 never fixed for now - let the solver calculate it from u_{0,1} ... u_{0, m} */
328
329 /* create CSC <-> COO exchange, init jacobians */
330 ✗ info.exc_jac = std::make_unique<ExchangeJacobians>(info);
331
332 /* create HESSIAN_PATTERNs and allocate buffers for extrapolation / evaluation */
333 ✗ info.exc_hes = std::make_unique<ExchangeHessians>(info);
334
335 /* create blocks (contains sparse patterns and mapping to buffer indices) */
336 ✗ GDOP::BoundarySweepLayout layout_mr(info.mayer_exists, info.r_size);
337 ✗ GDOP::FullSweepLayout layout_lfg(info.lagrange_exists, info.f_size, info.g_size);
338
339 /* fill layout_lfg and layout_mr objects with COO sparsity patterns */
340 ✗ init_eval(info, layout_lfg, layout_mr);
341 ✗ init_jac(info, layout_lfg, layout_mr);
342 ✗ init_hes(info, layout_lfg, layout_mr);
343
344 auto pc = std::make_unique<GDOP::ProblemConstants>(
345 ✗ info.mayer_exists,
346 ✗ info.lagrange_exists,
347 std::move(x_bounds),
348 std::move(u_bounds),
349 std::move(p_bounds),
350 std::move(T_bounds),
351 std::move(xu0_fixed),
352 std::move(xuf_fixed),
353 std::move(T_fixed),
354 std::move(r_bounds),
355 std::move(g_bounds),
356 mesh
357 ✗ );
358
359 ✗ auto fs = std::make_unique<FullSweep>(std::move(layout_lfg), *pc, info);
360 ✗ auto bs = std::make_unique<BoundarySweep>(std::move(layout_mr), *pc, info);
361 ✗ auto dyn = std::make_unique<Dynamics>(*pc, info);
362
363 ✗ return GDOP::Problem(std::move(fs), std::move(bs), std::move(pc), std::move(dyn));
364 ✗ }
365
366 } // namespace OpenModelica
367