Linux GNU 11.4.0 Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 78.7% 448 / 0 / 569
Functions: -% 0 / 1 / 1
Branches: 77.0% 208 / 0 / 270
OMCompiler/Compiler/BackEnd/ExpressionSolve.mo
Line Branch Exec Source
1 /*
2 * This file is part of OpenModelica.
3 *
4 * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC),
5 * c/o Linköpings universitet, Department of Computer and Information Science,
6 * SE-58183 Linköping, Sweden.
7 *
8 * All rights reserved.
9 *
10 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF AGPL VERSION 3 LICENSE OR
11 * THIS OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8.
12 * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES
13 * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GNU AGPL
14 * VERSION 3, ACCORDING TO RECIPIENTS CHOICE.
15 *
16 * The OpenModelica software and the OSMC (Open Source Modelica Consortium)
17 * Public License (OSMC-PL) are obtained from OSMC, either from the above
18 * address, from the URLs:
19 * http://www.openmodelica.org or
20 * https://github.com/OpenModelica/ or
21 * http://www.ida.liu.se/projects/OpenModelica,
22 * and in the OpenModelica distribution.
23 *
24 * GNU AGPL version 3 is obtained from:
25 * https://www.gnu.org/licenses/licenses.html#GPL
26 *
27 * This program is distributed WITHOUT ANY WARRANTY; without
28 * even the implied warranty of MERCHANTABILITY or FITNESS
29 * FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH
30 * IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF OSMC-PL.
31 *
32 * See the full OSMC Public License conditions for more details.
33 *
34 */
35
36 encapsulated package ExpressionSolve
37 " file: ExpressionSolve.mo
38 package: ExpressionSolve
39 description: ExpressionSolve
40
41
42 This file contains the module ExpressionSolve, which contains functions
43 to solve a DAE.Exp for a DAE.Exp"
44
45 // public imports
46 public import Absyn;
47 public import AbsynUtil;
48 public import DAE;
49
50 // protected imports
51 protected import ComponentReference;
52 protected import ComponentReferenceBasics;
53 protected import Debug;
54 protected import Differentiate;
55 protected import ElementSource;
56 protected import Ceval;
57 protected import Expression;
58 protected import ExpressionBasics;
59 protected import ExpressionDump;
60 protected import ExpressionSimplify;
61 protected import Flags;
62 protected import Global;
63 protected import List;
64 protected import Inline;
65 protected import BackendDAE;
66 protected import BackendDAEUtil;
67 protected import BackendEquation;
68 protected import BackendVariable;
69 protected import Types;
70
71 // =============================================================================
72 // section for postOptModule >>solveSimpleEquations<<
73 //
74 // solve simple equations otherwise detect EQUATIONSYSTEM
75 // =============================================================================
76
77 public function solveSimpleEquations
78 input output BackendDAE.BackendDAE dae;
79 protected
80 list<BackendDAE.EqSystem> systs;
81 BackendDAE.Shared shared;
82 algorithm
83 3936 (systs, shared) := List.mapFold(dae.eqs, solveSimpleEquationsSyst, dae.shared);
84 3936 dae := BackendDAE.DAE(systs, shared);
85 end solveSimpleEquations;
86
87 protected function solveSimpleEquationsSyst
88 input output BackendDAE.EqSystem syst;
89 input output BackendDAE.Shared shared;
90 protected
91 BackendDAE.StrongComponents comps = {}, oldComps;
92 BackendDAE.StrongComponent tmpComp;
93 array<Integer> ass1, ass2;
94 BackendDAE.Equation eqn;
95 BackendDAE.Var var;
96 Integer eindex, vindx;
97 Boolean solved;
98 algorithm
99 () := match syst
100 case BackendDAE.EQSYSTEM(matching = BackendDAE.MATCHING(ass1 = ass1, ass2 = ass2, comps = oldComps))
101 algorithm
102
2/2
✓ Branch 0 taken 160390 times.
✓ Branch 1 taken 53246 times.
 213636 for comp in oldComps loop
103 tmpComp := comp;
104
2/2
✓ Branch 1 taken 156026 times.
✓ Branch 2 taken 4364 times.
 160390 if isSingleEquation(comp) then
105
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 156026 times.
 156026 BackendDAE.SINGLEEQUATION(eqn=eindex, var=vindx) := comp;
106 156026 eqn := BackendEquation.get(syst.orderedEqs, eindex);
107
2/2
✓ Branch 1 taken 155777 times.
✓ Branch 2 taken 249 times.
 156026 if BackendEquation.isEquation(eqn) then
108 155777 var := BackendVariable.getVarAt(syst.orderedVars, vindx);
109 155777 (eqn, shared, solved) := solveSimpleEquation(eqn, var, shared);
110 155777 syst.orderedEqs := BackendEquation.setAtIndex(syst.orderedEqs, eindex, eqn);
111
2/2
✓ Branch 0 taken 358 times.
✓ Branch 1 taken 155419 times.
 155777 if not solved then
112 358 tmpComp := BackendDAE.EQUATIONSYSTEM({eindex}, {vindx}, BackendDAE.EMPTY_JACOBIAN(), BackendDAE.JAC_NONLINEAR(), false);
113 end if;
114 end if;
115 end if;
116 comps := tmpComp :: comps;
117 end for;
118 106492 syst.matching := BackendDAE.MATCHING(ass1, ass2, listReverse(comps));
119 then ();
120
121 else ();
122 end match;
123 end solveSimpleEquationsSyst;
124
125 protected function isSingleEquation
126 input BackendDAE.StrongComponent comp;
127 output Boolean b;
128 algorithm
129 b := match comp
130 case BackendDAE.SINGLEEQUATION() then true;
131 else false;
132 end match;
133 end isSingleEquation;
134
135 protected function solveSimpleEquation
136 input output BackendDAE.Equation eqn;
137 input BackendDAE.Var var "solve eqn with respect to var";
138 input output BackendDAE.Shared shared;
139 output Boolean solved;
140 protected
141 DAE.ComponentRef cr;
142 DAE.Exp lhs,rhs,e1,e2,varexp,e;
143 BackendDAE.EquationAttributes attr;
144 DAE.ElementSource source;
145 Boolean isContinuousIntegration = BackendDAEUtil.isSimulationDAE(shared);
146 algorithm
147
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 155777 times.
 155777 BackendDAE.EQUATION(exp=lhs, scalar=rhs, source=source, attr=attr) := eqn;
148 155777 (e1, e2) := (lhs, rhs);
149 155777 BackendDAE.VAR(varName = cr) := var;
150 155777 varexp := Expression.crefExp(cr);
151
2/2
✓ Branch 1 taken 3530 times.
✓ Branch 2 taken 152247 times.
 155777 if BackendVariable.isStateVar(var) then
152 3530 varexp := Expression.expDer(varexp);
153 3530 cr := ComponentReference.crefPrefixDer(cr);
154 end if;
155
156 // phi: aren't the types of e1 and e2 the same? Can we make the equation have a type?
157
3/4
✓ Branch 2 taken 151345 times.
✓ Branch 3 taken 4432 times.
✓ Branch 6 taken 151345 times.
✗ Branch 7 not taken.
 155777 if (Types.isIntegerOrRealOrSubTypeOfEither(Expression.typeof(e1)) and Types.isIntegerOrRealOrSubTypeOfEither(Expression.typeof(e2))) then
158 151345 (e1, e2) := preprocessingSolve(e1, e2, varexp, NONE(), SOME(shared.functionTree), NONE(), 0, false);
159 end if;
160
161 try
162 155777 e := solve2(e1, e2, varexp, SOME(shared.functionTree), NONE(), false, isContinuousIntegration);
163 155419 source := ElementSource.addSymbolicTransformationSolve(true, source, cr, e1, e2, e, {});
164 155419 eqn := BackendEquation.generateEquation(varexp, e, source, attr);
165 solved := true;
166 else
167 // only return new eqn if it can be solved explicitely because intermediate results can be numerically bad
168 // solves ticket #4293
169 // ToDo: do other preprocessing like multiplying by divisors?
170 solved := false;
171 end try;
172
173
4/4
✓ Branch 0 taken 358 times.
✓ Branch 1 taken 155419 times.
✓ Branch 3 taken 144158 times.
✓ Branch 4 taken 11261 times.
 155777 if solved and not isSolvedFor(lhs, rhs, varexp) then
174 11261 shared := checkSolveCoefficient(lhs, rhs, cr, source, shared);
175 end if;
176 end solveSimpleEquation;
177
178 protected function isSolvedFor
179 "Whether lhs = rhs already has the form x = f(..) for varexp, up to a sign."
180 input DAE.Exp lhs;
181 input DAE.Exp rhs;
182 input DAE.Exp varexp;
183 output Boolean b = true;
184 algorithm
185 try
186 155419 solveSimple(lhs, rhs, varexp, 0);
187 else
188 try
189 16019 solveSimple(rhs, lhs, varexp, 0);
190 else
191 b := false;
192 end try;
193 end try;
194 end isSolvedFor;
195
196 protected function checkSolveCoefficient
197 "Checks the coefficient that solving lhs = rhs for cr divides by. If it only
198 depends on parameters and is zero for their compile-time values, the
199 parameters are added to Global.structuralParameters, so that translateModel
200 evaluates them and translates the model again. Otherwise an assert that the
201 coefficient stays nonzero is added to the parameter asserts."
202 input DAE.Exp lhs;
203 input DAE.Exp rhs;
204 input DAE.ComponentRef cr "$DER-prefixed for a state";
205 input DAE.ElementSource source;
206 input output BackendDAE.Shared shared;
207 protected
208 DAE.Exp coef, value, cond;
209 list<DAE.ComponentRef> params, zeroParams;
210 list<String> names;
211 BackendDAE.Var v;
212 DAE.Type ty;
213 String msg;
214 DAE.Statement stmt;
215 algorithm
216 try
217 11261 coef := Differentiate.differentiateExpSolve(Expression.replaceDerOpInExp(Expression.expSub(lhs, rhs)), cr, SOME(shared.functionTree));
218 11259 (coef, _) := ExpressionSimplify.simplify(coef);
219
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 11259 times.
 11259 false := Types.isArray(Expression.typeof(coef));
220 11259 params := List.uniqueOnTrue(Expression.extractCrefsFromExp(coef), ComponentReferenceBasics.crefEqual);
221
2/2
✓ Branch 0 taken 7994 times.
✓ Branch 1 taken 3265 times.
 11259 false := listEmpty(params);
222
2/2
✓ Branch 0 taken 4165 times.
✓ Branch 1 taken 1767 times.
 5932 for p in params loop
223
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 3257 times.
 4165 (v :: _, _) := BackendVariable.getVar(p, shared.globalKnownVars);
224
4/4
✓ Branch 1 taken 2741 times.
✓ Branch 2 taken 516 times.
✓ Branch 4 taken 2667 times.
✓ Branch 5 taken 74 times.
 3257 true := BackendVariable.isParam(v) and BackendVariable.varFixed(v);
225 end for;
226 1767 value := evaluateParameterExp(coef, shared.globalKnownVars);
227
2/2
✓ Branch 1 taken 285 times.
✓ Branch 2 taken 1482 times.
 1767 true := Expression.isConst(value);
228 else
229 9779 return;
230 end try;
231
232
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1480 times.
 1482 if Expression.isZero(value) then
233
5/6
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2 times.
 4 zeroParams := list(p for p guard Expression.isZero(evaluateParameterExp(Expression.crefExp(p), shared.globalKnownVars)) in params);
234
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if listEmpty(zeroParams) then
235 zeroParams := params;
236 end if;
237 2 names := getGlobalRoot(Global.structuralParameters);
238
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 4 for p in zeroParams loop
239 2 names := List.unionElt(ComponentReferenceBasics.printComponentRefStr(ComponentReference.crefStripSubs(p)), names);
240 end for;
241 2 setGlobalRoot(Global.structuralParameters, names);
242 else
243 1480 ty := Expression.typeof(coef);
244 1480 cond := DAE.RELATION(coef, DAE.NEQUAL(ty), Expression.makeConstZero(ty), -1, NONE());
245
2/2
✓ Branch 1 taken 339 times.
✓ Branch 2 taken 1141 times.
 1480 msg := if Expression.isCref(coef)
246 then "Parameter " + ExpressionBasics.printExpStr(coef) + " cannot be set to zero at run time, because that causes a structural change in the equations; set " + ExpressionBasics.printExpStr(coef) + " = 0 in the model and recompile it."
247 else "The parameters in " + ExpressionBasics.printExpStr(coef) + " cannot be set so that it is zero at run time, because that causes a structural change in the equations; set their values in the model and recompile it.";
248 1480 stmt := DAE.STMT_ASSERT(cond, DAE.SCONST(msg), DAE.ASSERTIONLEVEL_ERROR, source);
249
2/2
✓ Branch 3 taken 256 times.
✓ Branch 4 taken 1224 times.
 1480 if not List.any(shared.parameterAsserts, function assertCondEqual(stmt2 = stmt)) then
250 2448 shared.parameterAsserts := stmt :: shared.parameterAsserts;
251 end if;
252 end if;
253 end checkSolveCoefficient;
254
255 protected function evaluateParameterExp
256 input DAE.Exp exp;
257 input BackendDAE.Variables globalKnownVars;
258 output DAE.Exp value;
259 algorithm
260 1769 (value, _) := Expression.traverseExpBottomUp(exp, BackendDAEUtil.replaceVarWithValue, globalKnownVars);
261 1769 (value, _) := ExpressionSimplify.simplify(value);
262 end evaluateParameterExp;
263
264 public function assertCondEqual
265 input DAE.Statement stmt1;
266 input DAE.Statement stmt2;
267 output Boolean b;
268 algorithm
269 b := match (stmt1, stmt2)
270 10103 case (DAE.STMT_ASSERT(), DAE.STMT_ASSERT()) then ExpressionBasics.expEqual(stmt1.cond, stmt2.cond);
271 else false;
272 end match;
273 end assertCondEqual;
274
275 protected function printTryToSolve
276 "for debugging"
277 input String instanceName "getInstanceName from caller";
278 input DAE.Exp inExp1 "lhs";
279 input DAE.Exp inExp2 "rhs";
280 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
281 algorithm
282 ✗ print(instanceName + " tries to solve: " +
283 ExpressionBasics.printExpStr(inExp1) + " = " + ExpressionBasics.printExpStr(inExp2) +
284 "\nwith respect to: " + ExpressionBasics.printExpStr(inExp3) + "\n");
285 end printTryToSolve;
286
287 public function solve
288 "Solves an equation consisting of a right hand side (rhs) and a
289 left hand side (lhs), with respect to the expression given as
290 third argument, usually a variable."
291 input DAE.Exp inExp1 "lhs";
292 input DAE.Exp inExp2 "rhs";
293 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
294 input Option<AvlTreePathFunction.Tree> functions = NONE() "need for solve modelica functions";
295 output DAE.Exp outExp;
296 output list<DAE.Statement> outAsserts;
297 protected
298 list<BackendDAE.Equation> dummy1;
299 list<DAE.ComponentRef> dummy2;
300 Integer dummyI;
301 algorithm
302 //printTryToSolve(getInstanceName(), inExp1, inExp2, inExp3);
303
304 (outExp,outAsserts,dummy1, dummy2, dummyI) := matchcontinue inExp1
305 92325 case _ then solveSimple(inExp1, inExp2, inExp3, 0);
306 11593 case _ then solveSimple(inExp2, inExp1, inExp3, 0);
307 10513 case _ then solveWork(inExp1, inExp2, inExp3, NONE(), functions, NONE(), 0, false, false);
308 else algorithm
309
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 3774 times.
 3775 if Flags.isSet(Flags.FAILTRACE) then
310 1 Error.addInternalError("Failed to solve \"" + ExpressionBasics.printExpStr(inExp1) + " = " + ExpressionBasics.printExpStr(inExp2) + "\" w.r.t. \"" + ExpressionBasics.printExpStr(inExp3) + "\"", sourceInfo());
311 end if;
312 3775 then fail();
313 end matchcontinue;
314
315 88550 (outExp,_) := ExpressionSimplify.simplify1(outExp);
316 end solve;
317
318
319 public function solve2
320 "Solves an equation with modelica function consisting of a right hand side (rhs) and a
321 left hand side (lhs), with respect to the expression given as
322 third argument, usually a variable.
323 "
324 input DAE.Exp inExp1 "lhs";
325 input DAE.Exp inExp2 "rhs";
326 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
327 input Option<AvlTreePathFunction.Tree> functions "need for solve modelica functions";
328 input Option<Integer> uniqueEqIndex "offset for tmp vars";
329 input Boolean doInline = true;
330 input Boolean isContinuousIntegration = false;
331 output DAE.Exp outExp;
332 output list<DAE.Statement> outAsserts;
333 output list<BackendDAE.Equation> eqnForNewVars "eqn for tmp vars";
334 output list<DAE.ComponentRef> newVarsCrefs;
335 protected
336 Integer dummyI;
337 algorithm
338 //printTryToSolve(getInstanceName(), inExp1, inExp2, inExp3);
339
340 (outExp,outAsserts,eqnForNewVars,newVarsCrefs,dummyI) := matchcontinue inExp1
341 340479 case _ then solveSimple(inExp1, inExp2, inExp3, 0);
342 13476 case _ then solveSimple(inExp2, inExp1, inExp3, 0);
343 12775 case _ then solveWork(inExp1, inExp2, inExp3, NONE(), functions, uniqueEqIndex, 0, doInline, isContinuousIntegration);
344 else algorithm
345
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 628 times.
 632 if Flags.isSet(Flags.FAILTRACE) then
346 4 Error.addInternalError("Failed to solve \"" + ExpressionBasics.printExpStr(inExp1) + " = " + ExpressionBasics.printExpStr(inExp2) + "\" w.r.t. \"" + ExpressionBasics.printExpStr(inExp3) + "\"", sourceInfo());
347 end if;
348 632 then fail();
349 end matchcontinue;
350 end solve2;
351
352
353 protected function solveWork
354 input DAE.Exp inExp1 "lhs";
355 input DAE.Exp inExp2 "rhs";
356 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
357 input Option<DAE.Exp> optCond "condition from an if expression";
358 input Option<AvlTreePathFunction.Tree> functions;
359 input Option<Integer> uniqueEqIndex "offset for tmp vars";
360 input Integer idepth;
361 input Boolean doInline;
362 input Boolean isContinuousIntegration;
363 output DAE.Exp outExp;
364 output list<DAE.Statement> outAsserts;
365 output list<BackendDAE.Equation> eqnForNewVars "eqn for tmp vars";
366 output list<DAE.ComponentRef> newVarsCrefs;
367 output Integer depth;
368 protected
369 DAE.Exp e1, e2;
370 list<BackendDAE.Equation> eqnForNewVars1, eqnForNewVars2;
371 list<DAE.ComponentRef> newVarsCrefs1, newVarsCrefs2;
372 algorithm
373 (e1, e2, eqnForNewVars1, newVarsCrefs1, depth) := matchcontinue inExp1
374 23932 case _ then preprocessingSolve(inExp1, inExp2, inExp3, optCond, functions, uniqueEqIndex, idepth, doInline);
375 else
376 algorithm
377
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 21 times.
 21 if Flags.isSet(Flags.FAILTRACE) then
378 ✗ Debug.trace("\n-ExpressionSolve.preprocessingSolve failed:\n");
379 ✗ Debug.trace(ExpressionBasics.printExpStr(inExp1) + " = " + ExpressionBasics.printExpStr(inExp2));
380 ✗ Debug.trace(" with respect to: " + ExpressionBasics.printExpStr(inExp3));
381 end if;
382 21 then (inExp1,inExp2,{},{}, idepth);
383 end matchcontinue;
384
385 (outExp, outAsserts, eqnForNewVars2, newVarsCrefs2, depth) := matchcontinue e1
386 23932 case _ then solveIfExp(e1, e2, inExp3, optCond, functions, uniqueEqIndex, depth, doInline, isContinuousIntegration);
387 23685 case _ then solveSimple(e1, e2, inExp3, depth);
388 4606 case _ then solveLinearSystem(e1, e2, inExp3, functions, depth);
389 end matchcontinue;
390
391 19427 eqnForNewVars := listAppend(eqnForNewVars1, eqnForNewVars2);
392 19427 newVarsCrefs := listAppend(newVarsCrefs1, newVarsCrefs2);
393 end solveWork;
394
395 protected function solveSimple
396 "Solves simple equations like
397 a = f(..)
398 der(a) = f(..)
399 -a = f(..)
400 -der(a) = f(..)"
401 input DAE.Exp inExp1 "lhs";
402 input DAE.Exp inExp2 "rhs";
403 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
404 input Integer idepth;
405 output DAE.Exp outExp;
406 output list<DAE.Statement> outAsserts;
407 output list<BackendDAE.Equation> eqnForNewVars = {} "eqn for tmp vars";
408 output list<DAE.ComponentRef> newVarsCrefs = {};
409 output Integer odepth = idepth;
410
411 algorithm
412 //printTryToSolve(getInstanceName(), inExp1, inExp2, inExp3);
413
414 (outExp,outAsserts) := match (inExp1, inExp3)
415 local
416 DAE.ComponentRef cr,cr1;
417 DAE.Type tp;
418 list<DAE.Statement> asserts;
419
420 // special case when already solved, cr1 = rhs, otherwise division by zero when dividing with derivative
421 case (DAE.CREF(componentRef = cr1), DAE.CREF(componentRef = cr))
422 guard ComponentReferenceBasics.crefEqual(cr, cr1) and (not Expression.expHasCrefNoPreOrStart(inExp2, cr))
423 then
424 (inExp2,{});
425 case (DAE.CALL(path = Absyn.IDENT(name = "der"),expLst = {DAE.CREF(componentRef = cr1)}), DAE.CALL(path = Absyn.IDENT(name = "der"),expLst = {DAE.CREF(componentRef = cr)}))
426 guard ComponentReferenceBasics.crefEqual(cr, cr1) and (not Expression.expHasDerCref(inExp2, cr))
427 then
428 (inExp2,{});
429
430 // -cr = exp
431 case (DAE.UNARY(operator = DAE.UMINUS(), exp = DAE.CREF(componentRef = cr1)), DAE.CREF(componentRef = cr))
432 guard ComponentReferenceBasics.crefEqual(cr1,cr) and (not Expression.expHasCrefNoPreOrStart(inExp2,cr))
433 4523 then
434 (Expression.negate(inExp2),{});
435 case (DAE.UNARY(operator = DAE.UMINUS_ARR(), exp = DAE.CREF(componentRef = cr1)), DAE.CREF(componentRef = cr))
436 guard ComponentReferenceBasics.crefEqual(cr1,cr) and (not Expression.expHasCrefNoPreOrStart(inExp2,cr)) // cr not in e2
437 ✗ then
438 (Expression.negate(inExp2),{});
439 case (DAE.UNARY(operator = DAE.UMINUS(), exp = DAE.CALL(path = Absyn.IDENT(name = "der"),expLst = {DAE.CREF(componentRef = cr1)})), DAE.CALL(path = Absyn.IDENT(name = "der"),expLst = {DAE.CREF(componentRef = cr)}))
440 guard ComponentReferenceBasics.crefEqual(cr1,cr) and (not Expression.expHasDerCref(inExp2,cr)) // cr not in e2
441 2 then
442 (Expression.negate(inExp2),{});
443 case (DAE.UNARY(operator = DAE.UMINUS_ARR(), exp = DAE.CALL(path = Absyn.IDENT(name = "der"),expLst = {DAE.CREF(componentRef = cr1)})), DAE.CALL(path = Absyn.IDENT(name = "der"),expLst = {DAE.CREF(componentRef = cr)}))
444 guard ComponentReferenceBasics.crefEqual(cr1,cr) and (not Expression.expHasDerCref(inExp2,cr))
445 ✗ then
446 (Expression.negate(inExp2),{});
447
448 // !cr = exp
449 case (DAE.LUNARY(operator = DAE.NOT(), exp = DAE.CREF(componentRef = cr1)), DAE.CREF(componentRef = cr))
450 guard ComponentReferenceBasics.crefEqual(cr1,cr) and (not Expression.expHasCrefNoPreOrStart(inExp2,cr))
451 266 then
452 (Expression.negate(inExp2),{});
453
454 // Integer(enumcr) = ...
455 case (DAE.CALL(path = Absyn.IDENT(name = "Integer"),expLst={DAE.CREF(componentRef = cr1)}), DAE.CREF(componentRef = cr,ty=tp))
456 guard ComponentReferenceBasics.crefEqual(cr, cr1) and (not Expression.expHasCrefNoPreorDer(inExp2,cr))
457 algorithm
458 ✗ asserts := generateAssertType(tp,cr,inExp3,{});
459 ✗ then (DAE.CAST(tp,inExp2),asserts);
460 else fail();
461 end match;
462 end solveSimple;
463
464 protected function generateAssertType
465 input DAE.Type tp;
466 input DAE.ComponentRef cr;
467 input DAE.Exp iExp;
468 input list<DAE.Statement> inAsserts;
469 output list<DAE.Statement> outAsserts;
470 algorithm
471 outAsserts := match tp
472 local
473 Absyn.Path path,p1,pn;
474 list<String> names;
475 Integer n;
476 DAE.Exp e1,en,e,es;
477 String s1,sn,estr,crstr;
478 case DAE.T_ENUMERATION(path=path,names=names)
479 algorithm
480 ✗ p1 := AbsynUtil.suffixPath(path,listHead(names));
481 ✗ e1 := DAE.ENUM_LITERAL(p1,1);
482 ✗ n := listLength(names);
483 ✗ pn := AbsynUtil.suffixPath(path,listGet(names,n));
484 ✗ en := DAE.ENUM_LITERAL(p1,n);
485 ✗ s1 := AbsynUtil.pathString(p1);
486 ✗ sn := AbsynUtil.pathString(pn);
487 ✗ crstr := ComponentReferenceBasics.printComponentRefStr(cr);
488 ✗ estr := "Expression for " + crstr + " out of min(" + s1 + ")/max(" + sn + ") = ";
489 // iExp >= e1 and iExp <= en
490 ✗ e := DAE.LBINARY(DAE.RELATION(iExp,DAE.GREATEREQ(DAE.T_INTEGER_DEFAULT),e1,-1,NONE()),DAE.AND(DAE.T_BOOL_DEFAULT),
491 DAE.RELATION(iExp,DAE.LESSEQ(DAE.T_INTEGER_DEFAULT),en,-1,NONE()));
492 ✗ es := Expression.makePureBuiltinCall("String", {iExp,DAE.SCONST("d")}, DAE.T_STRING_DEFAULT);
493 ✗ es := DAE.BINARY(DAE.SCONST(estr),DAE.ADD(DAE.T_STRING_DEFAULT),es);
494 ✗ then
495 DAE.STMT_ASSERT(e,es,DAE.ASSERTIONLEVEL_ERROR,DAE.emptyElementSource)::inAsserts;
496 else inAsserts;
497 end match;
498 end generateAssertType;
499
500 public function preprocessingSolve
501 "
502 preprocessing for solve1,
503 sorting and split terms , with respect to the expression given as
504 third argument.
505
506 {f(x,y), g(x,y),x} -> {h(x), k(y)}
507
508 author: Vitalij Ruge
509 "
510
511 input output DAE.Exp x "lhs";
512 input output DAE.Exp y "rhs";
513 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
514 input Option<DAE.Exp> optCond "condition from an if expression";
515 input Option<AvlTreePathFunction.Tree> functions;
516 input Option<Integer> uniqueEqIndex "offset for tmp vars";
517 input Integer idepth;
518 input Boolean doInline;
519 output list<BackendDAE.Equation> eqnForNewVars = {} "eqn for tmp vars";
520 output list<DAE.ComponentRef> newVarsCrefs = {};
521 output Integer depth = idepth;
522
523 protected
524 DAE.Exp lhsX, rhsX, lhsY, rhsY, N;
525 Boolean con, new_x, inlineFun = true;
526 Integer iter;
527 Integer numSimplifed = 0 ;
528 algorithm
529 // split and sort
530 175277 (lhsX, lhsY) := preprocessingSolve5(x, inExp3,true);
531 175277 (rhsX, rhsY) := preprocessingSolve5(y, inExp3,true);
532 175277 x := Expression.expSub(lhsX, rhsX);
533 175277 y := Expression.expSub(rhsY, lhsY);
534
535 175256 con := not Expression.isCref(x);
536 iter := 0;
537
2/2
✓ Branch 0 taken 35531 times.
✓ Branch 1 taken 139725 times.
 175256 if con then
538 35531 x := unifyFunCalls(x, inExp3);
539 end if;
540
4/4
✓ Branch 0 taken 41373 times.
✓ Branch 1 taken 148693 times.
✓ Branch 3 taken 40759 times.
✓ Branch 4 taken 614 times.
 190066 while con and iter < 1000 and (not Expression.isCref(x)) loop
541 40759 (x, y, con) := preprocessingSolve2(x,y, inExp3);
542 40759 (x, y, new_x) := preprocessingSolve3(x,y, inExp3);
543
4/4
✓ Branch 0 taken 10154 times.
✓ Branch 1 taken 30605 times.
✓ Branch 2 taken 10121 times.
✓ Branch 3 taken 33 times.
 50880 con := con or new_x;
544
2/2
✓ Branch 0 taken 61 times.
✓ Branch 1 taken 40759 times.
 40820 while new_x loop
545 61 (x, y, new_x) := preprocessingSolve3(x,y, inExp3);
546 end while;
547
548
2/2
✓ Branch 1 taken 14810 times.
✓ Branch 2 taken 25949 times.
 40759 if Expression.isCref(x) then
549 break;
550 end if;
551 14810 (x, y, new_x) := removeSimpleCalls(x,y, inExp3);
552
4/4
✓ Branch 0 taken 10121 times.
✓ Branch 1 taken 4689 times.
✓ Branch 2 taken 10070 times.
✓ Branch 3 taken 51 times.
 14810 con := con or new_x;
553 14810 (x, y, new_x) := preprocessingSolve4(x,y, inExp3);
554
4/4
✓ Branch 0 taken 14803 times.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 10063 times.
✓ Branch 3 taken 4740 times.
 14810 con := new_x or con;
555 // TODO: use new defined function, which missing in the cpp runtime
556
7/8
✗ Branch 0 not taken.
✓ Branch 1 taken 14810 times.
✓ Branch 2 taken 1821 times.
✓ Branch 3 taken 12989 times.
✓ Branch 5 taken 70 times.
✓ Branch 6 taken 1751 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 68 times.
 14810 if isSome(uniqueEqIndex) and not stringEqual(Config.simCodeTarget(), "Cpp") then
557 1753 (x, y, new_x, eqnForNewVars, newVarsCrefs, depth) := preprocessingSolveTmpVars(x, y, inExp3, optCond, Util.getOption(uniqueEqIndex), eqnForNewVars, newVarsCrefs, depth);
558
4/4
✓ Branch 0 taken 1494 times.
✓ Branch 1 taken 259 times.
✓ Branch 2 taken 962 times.
✓ Branch 3 taken 532 times.
 2715 con := new_x or con;
559 end if;
560
561
2/2
✓ Branch 0 taken 9887 times.
✓ Branch 1 taken 4923 times.
 14810 if (not con) then
562
2/2
✓ Branch 0 taken 9821 times.
✓ Branch 1 taken 66 times.
 9887 if (numSimplifed < 3) then
563 9821 (x, con) := ExpressionSimplify.simplify(x);
564 9821 numSimplifed := numSimplifed + 1;
565 end if;
566 // Z/N = rhs -> Z = rhs*N
567 9887 (x,N) := Expression.makeFraction(x);
568
2/2
✓ Branch 1 taken 279 times.
✓ Branch 2 taken 9608 times.
 9887 if not Expression.isOne(N) then
569 //print("\nx ");print(ExpressionBasics.printExpStr(x));print("\nN ");print(ExpressionBasics.printExpStr(N));
570 279 new_x := true;
571 279 y := Expression.expMul(y,N);
572 end if;
573
574
4/4
✓ Branch 0 taken 9608 times.
✓ Branch 1 taken 279 times.
✓ Branch 2 taken 8999 times.
✓ Branch 3 taken 609 times.
 18886 con := new_x or con;
575 end if;
576
577
2/2
✓ Branch 0 taken 5811 times.
✓ Branch 1 taken 8999 times.
 14810 if con then
578 5811 (lhsX, lhsY) := preprocessingSolve5(x, inExp3, true);
579 5811 (rhsX, rhsY) := preprocessingSolve5(y, inExp3, false);
580 5811 x := Expression.expSub(lhsX, rhsX);
581 5811 y := Expression.expSub(rhsY, lhsY);
582 elseif doInline and inlineFun then
583 507 iter := iter + 50;
584 if inlineFun then
585 507 (x,con) := solveFunCalls(x, inExp3, functions);
586 inlineFun := false;
587
2/2
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 476 times.
 507 if con then
588 numSimplifed := 0;
589 end if;
590 end if;
591 end if;
592
593 14810 iter := iter + 1;
594 //print("\nx ");print(ExpressionBasics.printExpStr(x));print("\ny ");print(ExpressionBasics.printExpStr(y));
595 end while;
596
597 175256 y := ExpressionSimplify.simplify1(y);
598 end preprocessingSolve;
599
600 protected function preprocessingSolve2
601 "
602 helprer function for preprocessingSolve
603 e.g.
604 x/(x+c1) = -c2 --> x + (x+c1)*c2 = 0
605
606 author: Vitalij Ruge
607 "
608 input DAE.Exp inExp1 "lhs";
609 input DAE.Exp inExp2 "rhs";
610 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
611
612 output DAE.Exp olhs;
613 output DAE.Exp orhs;
614 output Boolean con "continue";
615
616 algorithm
617
618 (olhs, orhs, con) := match inExp1
619 local
620 DAE.Exp e,b, fa, ga, lhs;
621 DAE.Type tp;
622 list<DAE.Exp> eWithX, factorWithX, factorWithoutX;
623 DAE.Exp pWithX, pWithoutX;
624
625 // -f(a) = b => f(a) = -b
626 case DAE.UNARY(DAE.UMINUS(), fa)
627 guard expHasCref(fa, inExp3) and not expHasCref(inExp2, inExp3)
628 algorithm
629 18529 b := Expression.negate(inExp2);
630 then (fa, b, true);
631
632 case DAE.UNARY(DAE.UMINUS_ARR(), fa)
633 guard expHasCref(fa, inExp3) and not expHasCref(inExp2, inExp3)
634 algorithm
635 ✗ b := Expression.negate(inExp2);
636 then (fa, b, true);
637
638 // b/f(a) = rhs => f(a) = b/rhs solve for a
639 case DAE.BINARY(b,DAE.DIV(_),fa)
640 guard expHasCref(fa, inExp3) and (not expHasCref(b, inExp3)) and (not expHasCref(inExp2, inExp3))
641 algorithm
642 31 e := Expression.makeDiv(b, inExp2);
643 then(fa, e, true);
644
645 // b*f(a) = rhs => f(a) = rhs/b solve for a
646 case DAE.BINARY(b, DAE.MUL(_), fa)
647 guard expHasCref(fa, inExp3) and (not expHasCref(b, inExp3)) and (not expHasCref(inExp2, inExp3))
648 algorithm
649
650 11166 eWithX := Expression.expandFactors(inExp1);
651 11166 (factorWithX, factorWithoutX) := List.split1OnTrue(eWithX, expHasCref, inExp3);
652 11166 pWithX := makeProductLstSort(factorWithX);
653 11166 pWithoutX := makeProductLstSort(factorWithoutX);
654
655 11166 e := Expression.makeDiv(inExp2, pWithoutX);
656
657 then(pWithX, e, true);
658
659 // b*a = rhs => a = rhs/b solve for a
660 case DAE.BINARY(b, DAE.MUL(_), fa)
661 guard expHasCref(fa, inExp3) and (not expHasCref(b, inExp3)) and (not expHasCref(inExp2, inExp3))
662 algorithm
663 ✗ e := Expression.makeDiv(inExp2, b);
664 then(fa, e, true);
665
666 // a*b = rhs => a = rhs/b solve for a
667 case DAE.BINARY(fa, DAE.MUL(_), b)
668 guard expHasCref(fa, inExp3) and (not expHasCref(b, inExp3)) and (not expHasCref(inExp2, inExp3))
669 algorithm
670 ✗ e := Expression.makeDiv(inExp2, b);
671 then(fa, e, true);
672
673 // f(a)/b = rhs => f(a) = rhs*b solve for a
674 case DAE.BINARY(fa, DAE.DIV(_), b)
675 guard expHasCref(fa, inExp3) and (not expHasCref(b, inExp3)) and (not expHasCref(inExp2, inExp3))
676 algorithm
677 869 e := Expression.expMul(inExp2, b);
678 then (fa, e, true);
679
680 // g(a)/f(a) = rhs => rhs*f(a) - g(a) = 0 solve for a
681 case DAE.BINARY(ga, DAE.DIV(tp), fa)
682 guard expHasCref(fa, inExp3) and expHasCref(ga, inExp3) and (not expHasCref(inExp2, inExp3))
683 algorithm
684
685 10 e := Expression.expMul(inExp2, fa);
686 10 lhs := Expression.expSub(e, ga);
687 10 e := Expression.makeConstZero(tp);
688
689 then(lhs, e, true);
690
691 else (inExp1, inExp2, false);
692
693 end match;
694
695 end preprocessingSolve2;
696
697 protected function preprocessingSolve3
698 "
699 helprer function for preprocessingSolve
700
701 (r1)^f(a) = r2 => f(a) = ln(r2)/ln(r1)
702 f(a)^b = 0 => f(a) = 0
703 f(a)^n = c => f(a) = c^(1/n)
704 abs(x) = 0
705 author: Vitalij Ruge
706 "
707 input DAE.Exp inExp1 "lhs";
708 input DAE.Exp inExp2 "rhs";
709 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
710
711 output DAE.Exp olhs;
712 output DAE.Exp orhs;
713 output Boolean con "continue";
714
715 algorithm
716 (olhs, orhs, con) := match(inExp1, inExp2)
717 local
718 Real r, r1, r2;
719 DAE.Exp e1, e2, res;
720
721 // (r1)^f(a) = r2 => f(a) = ln(r2)/ln(r1)
722 case (DAE.BINARY(e1 as DAE.RCONST(r1),DAE.POW(_),e2), DAE.RCONST(r2))
723 guard r2 > 0.0 and r1 > 0.0 and (not Expression.isConstOne(e1)) and expHasCref(e2, inExp3)
724 algorithm
725 ✗ r := log(r2) / log(r1);
726 ✗ res := DAE.RCONST(r);
727 then
728 (e2, res, true);
729
730 // f(a)^b = 0 => f(a) = 0
731 case (DAE.BINARY(e1,DAE.POW(_),e2), DAE.RCONST(real = 0.0))
732 guard expHasCref(e1, inExp3) and (not expHasCref(e2, inExp3))
733 then
734 (e1, inExp2, true);
735
736 // f(a)^n = c => f(a) = c^(1/n)
737 // where n is odd
738 case (DAE.BINARY(e1,DAE.POW(_),e2 as DAE.RCONST(r)), _)
739 guard (not expHasCref(inExp2, inExp3)) and expHasCref(e1, inExp3) and (1.0 == realMod(r,2.0))
740 algorithm
741 40 res := Expression.makeDiv(DAE.RCONST(1.0),e2);
742 40 res := Expression.expPow(inExp2,res);
743 then
744 (e1, res, true);
745
746 // sqrt(f(a)) = f(a)^n = c => f(a) = c^(1/n)
747 case (DAE.BINARY(e1,DAE.POW(_),DAE.RCONST(0.5)), _)
748 guard not expHasCref(inExp2, inExp3) and expHasCref(e1, inExp3)
749 algorithm
750 4 res := Expression.expPow(inExp2,DAE.RCONST(2.0));
751 then
752 (e1, res, true);
753
754 // abs(x) = 0
755 case (DAE.CALL(path = Absyn.IDENT(name = "abs"),expLst = {e1}), DAE.RCONST(0.0))
756 then (e1,inExp2,true);
757
758 // sign(x) = 0
759 case (DAE.CALL(path = Absyn.IDENT(name = "sign"),expLst = {e1}), DAE.RCONST(0.0))
760 then (e1,inExp2,true);
761
762
763 else (inExp1, inExp2, false);
764
765 end match;
766
767
768 end preprocessingSolve3;
769
770
771 protected function preprocessingSolve4
772
773 "
774 helprer function for preprocessingSolve
775
776 e.g.
777 sqrt(f(x)) - sqrt(g(x))) = 0 = f(x) - g(x)
778 exp(f(x)) - exp(g(x))) = 0 = f(x) - g(x)
779
780 author: Vitalij Ruge
781 "
782
783 input DAE.Exp inExp1;
784 input DAE.Exp inExp2;
785 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
786 output DAE.Exp oExp1;
787 output DAE.Exp oExp2;
788 output Boolean newX;
789
790 algorithm
791
792 (oExp1, oExp2, newX) := match(inExp1, inExp2)
793 local
794 DAE.Exp e1,e2,e3,e4, e, e_1, e_2;
795 DAE.Type tp;
796
797 // exp(f(x)) - exp(g(x)) = 0
798 case(DAE.BINARY(DAE.CALL(path = Absyn.IDENT("exp"), expLst={e1}), DAE.SUB(_),
799 DAE.CALL(path = Absyn.IDENT("exp"), expLst={e2})), DAE.RCONST(0.0))
800 then (e1, e2, true);
801
802 // log(f(x)) - log(g(x)) = 0
803 case(DAE.BINARY(DAE.CALL(path = Absyn.IDENT("log"), expLst={e1}), DAE.SUB(_),
804 DAE.CALL(path = Absyn.IDENT("log"), expLst={e2})), DAE.RCONST(0.0))
805 then (e1, e2, true);
806
807 // log10(f(x)) - log10(g(x)) = 0
808 case(DAE.BINARY(DAE.CALL(path = Absyn.IDENT("log10"), expLst={e1}), DAE.SUB(_),
809 DAE.CALL(path = Absyn.IDENT("log10"), expLst={e2})), DAE.RCONST(0.0))
810 then (e1, e2, true);
811
812 // sinh(f(x)) - sinh(g(x)) = 0
813 case(DAE.BINARY(DAE.CALL(path = Absyn.IDENT("sinh"), expLst={e1}), DAE.SUB(_),
814 DAE.CALL(path = Absyn.IDENT("sinh"), expLst={e2})), DAE.RCONST(0.0))
815 then (e1, e2, true);
816
817 // tanh(f(x)) - tanh(g(x)) = 0
818 case(DAE.BINARY(DAE.CALL(path = Absyn.IDENT("tanh"), expLst={e1}), DAE.SUB(_),
819 DAE.CALL(path = Absyn.IDENT("tanh"), expLst={e2})), DAE.RCONST(0.0))
820 then (e1, e2, true);
821
822 // sqrt(f(x)) - sqrt(g(x)) = 0
823 case(DAE.BINARY(DAE.CALL(path = Absyn.IDENT("sqrt"), expLst={e1}), DAE.SUB(_),
824 DAE.CALL(path = Absyn.IDENT("sqrt"), expLst={e2})), DAE.RCONST(0.0))
825 then (e1, e2, true);
826
827 // sinh(f(x)) - cosh(g(x)) = 0
828 case(DAE.BINARY(DAE.CALL(path = Absyn.IDENT("sinh"), expLst={e1}), DAE.SUB(_),
829 DAE.CALL(path = Absyn.IDENT("cosh"), expLst={e2})), DAE.RCONST(0.0))
830 guard ExpressionBasics.expEqual(e1,e2)
831 then (e1, inExp2, true);
832
833 case(DAE.BINARY(DAE.CALL(path = Absyn.IDENT("cosh"), expLst={e1}), DAE.SUB(_),
834 DAE.CALL(path = Absyn.IDENT("sinh"), expLst={e2})), DAE.RCONST(0.0))
835 guard ExpressionBasics.expEqual(e1,e2)
836 then (e1, inExp2, true);
837
838
839 // y*sinh(x) - z*cosh(x) = 0
840 case(DAE.BINARY(DAE.BINARY(e3,DAE.MUL(),DAE.CALL(path = Absyn.IDENT("sinh"), expLst={e1})), DAE.SUB(tp),
841 DAE.BINARY(e4,DAE.MUL(),DAE.CALL(path = Absyn.IDENT("cosh"), expLst={e2}))), DAE.RCONST(0.0))
842 guard ExpressionBasics.expEqual(e1,e2)
843 algorithm
844 2 e := Expression.makePureBuiltinCall("tanh",{e1},tp);
845 2 then (Expression.expMul(e3,e), e4, true);
846
847 case(DAE.BINARY(DAE.BINARY(e4,DAE.MUL(),DAE.CALL(path = Absyn.IDENT("cosh"), expLst={e2})), DAE.SUB(tp),
848 DAE.BINARY(e3,DAE.MUL(),DAE.CALL(path = Absyn.IDENT("sinh"), expLst={e1}))), DAE.RCONST(0.0))
849 guard ExpressionBasics.expEqual(e1,e2)
850 algorithm
851 ✗ e := Expression.makePureBuiltinCall("tanh",{e1},tp);
852 ✗ then (Expression.expMul(e3,e), e4, true);
853
854
855
856 // sqrt(x) - x = 0 -> x = x^2
857 case(DAE.BINARY(DAE.CALL(path = Absyn.IDENT("sqrt"), expLst={e1}), DAE.SUB(_),e2), DAE.RCONST(0.0))
858 ✗ then (e1, Expression.expPow(e2, DAE.RCONST(2.0)), true);
859
860 case(DAE.BINARY(e2, DAE.SUB(_),DAE.CALL(path = Absyn.IDENT("sqrt"), expLst={e1})), DAE.RCONST(0.0))
861 1 then (e1, Expression.expPow(e2, DAE.RCONST(2.0)), true);
862
863 // f(x)^n - g(x)^n = 0 -> (f(x)/g(x))^n = 1
864 case(DAE.BINARY(DAE.BINARY(e1, DAE.POW(), e2), DAE.SUB(tp), DAE.BINARY(e3, DAE.POW(), e4)), DAE.RCONST(0.0))
865 guard ExpressionBasics.expEqual(e2,e4) and expHasCref(e1,inExp3) and expHasCref(e3,inExp3)
866 algorithm
867 2 e := Expression.expPow(Expression.makeDiv(e1,e3),e2);
868 2 (e_1, e_2, _) := preprocessingSolve3(e, Expression.makeConstOne(tp), inExp3);
869 2 then (e_1, e_2, true);
870
871 else (inExp1, inExp2, false);
872
873 end match;
874
875
876 end preprocessingSolve4;
877
878 protected function expAddX
879 "
880 helprer function for preprocessingSolve
881
882 if(y,g(x),h(x)) + x => if(y, g(x) + x, h(x) + x)
883 a*f(x) + b*f(x) = (a+b)*f(x)
884 author: Vitalij Ruge
885 "
886 input DAE.Exp inExp1 "lhs";
887 input DAE.Exp inExp2 "rhs";
888 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
889
890 output DAE.Exp ores;
891
892 algorithm
893 ores := matchcontinue(inExp1, inExp2)
894 local
895 DAE.Exp e, e1, e2, e3, e4, res;
896
897 case(DAE.IFEXP(e,e1,e2), _)
898 guard expHasCref(e1, inExp3) and expHasCref(e2, inExp3) and (not expHasCref(e, inExp3))
899 algorithm
900 1036 e3 := expAddX(inExp2, e1, inExp3);
901 1036 e4 := expAddX(inExp2, e2, inExp3);
902
903 1036 res := DAE.IFEXP(e, e3, e4);
904 then res;
905
906 case(_, DAE.IFEXP(e,e1,e2))
907 guard expHasCref(e1, inExp3) and expHasCref(e2, inExp3) and (not expHasCref(e, inExp3))
908 algorithm
909 476 e3 := expAddX(inExp1, e1, inExp3);
910 476 e4 := expAddX(inExp1, e2, inExp3);
911
912 476 res := DAE.IFEXP(e, e3, e4);
913 then res;
914
915 else
916 algorithm
917 365202 res := expAddX2(inExp1, inExp2, inExp3);
918 then res;
919
920 end matchcontinue;
921
922 end expAddX;
923
924 protected function expAddX2
925 "
926 helprer function for preprocessingSolve
927 a*f(x) + b*f(x) = (a+b)*f(x)
928 author: Vitalij Ruge
929 "
930 input DAE.Exp inExp1 "lhs";
931 input DAE.Exp inExp2 "rhs";
932 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
933
934 output DAE.Exp ores;
935
936 protected
937 list<DAE.Exp> f1, f2;
938 DAE.Exp e0,e1,e2;
939 Boolean neg;
940 list<DAE.Exp> factorWithX1, factorWithoutX1, factorWithX2, factorWithoutX2;
941 DAE.Exp pWithX1, pWithoutX1, pWithX2, pWithoutX2;
942
943 algorithm
944 (e0, e1, neg) := match inExp1
945 local DAE.Exp ee1, ee2;
946 case DAE.BINARY(ee1,DAE.ADD(),ee2)
947 then(ee1, ee2, false);
948 case DAE.BINARY(ee1,DAE.SUB(),ee2)
949 then(ee1, ee2, true);
950 else
951 then(DAE.RCONST(0.0), inExp1, false);
952 end match;
953
954 365202 f1 := Expression.expandFactors(e1);
955 365202 (factorWithX1, factorWithoutX1) := List.split1OnTrue(f1, expHasCref, inExp3);
956 365202 pWithX1 := makeProductLstSort(factorWithX1);
957 365202 pWithoutX1 := makeProductLstSort(factorWithoutX1);
958 365202 f2 := Expression.expandFactors(inExp2);
959 365202 (factorWithX2, factorWithoutX2) := List.split1OnTrue(f2, expHasCref, inExp3);
960 365202 (pWithX2,_) := ExpressionSimplify.simplify1(makeProductLstSort(factorWithX2));
961 365202 pWithoutX2 := makeProductLstSort(factorWithoutX2);
962 //print("\nf1 =");print(ExpressionDump.printExpListStr(f1));
963 //print("\nf2 =");print(ExpressionDump.printExpListStr(f2));
964
965
2/2
✓ Branch 1 taken 3530 times.
✓ Branch 2 taken 361672 times.
 365202 if ExpressionBasics.expEqual(pWithX2,pWithX1) then
966 // e0 + a*x + b*x -> e0 + (a+b)*x
967
2/2
✓ Branch 0 taken 3525 times.
✓ Branch 1 taken 5 times.
 3530 if not neg then
968 3525 ores := Expression.expAdd(pWithoutX1, pWithoutX2);
969 else
970 // e0 - a*x + b*x -> e0 + (b-a)*x
971 5 ores := Expression.expSub(pWithoutX2, pWithoutX1);
972 end if;
973 3530 ores := Expression.expMul(ores, pWithX2);
974 elseif ExpressionBasics.expEqual(pWithX2, Expression.negate(pWithX1)) then
975 // e0 + a*(-x) + b*x -> e0 + (b-a)*x
976 ✗ if not neg then
977 ✗ ores := Expression.expSub(pWithoutX2, pWithoutX1);
978 else
979 // e0 - a*(-x) + b*x -> e0 + (b-a)*x
980 ✗ ores := Expression.expAdd(pWithoutX1, pWithoutX2);
981 end if;
982 ✗ ores := Expression.expMul(ores, pWithX2);
983 else
984 361672 e1 := Expression.expMul(pWithoutX1, pWithX1);
985 361672 e2 := Expression.expMul(pWithoutX2, pWithX2);
986 361672 ores := Expression.expAdd(e1,e2);
987 end if;
988
989 365202 ores := Expression.expAdd(e0,ores);
990
991 end expAddX2;
992
993 public function collectX
994 input DAE.Exp inExp1 "lhs";
995 input DAE.Exp inExp3 "DAE.CREF";
996 input Boolean expand = true;
997 output DAE.Exp outLhs;
998 output DAE.Exp outRhs;
999 algorithm
1000 15 (outLhs, outRhs) := preprocessingSolve5(inExp1, inExp3, expand);
1001 end collectX;
1002
1003 protected function preprocessingSolve5
1004 "
1005 helprer function for preprocessingSolve
1006 split and sort with respect to x
1007 where x = cref
1008
1009 f(x,y) = {h(y)*g(x,y), k(y)}
1010
1011 author: Vitalij Ruge
1012 "
1013 input DAE.Exp inExp1 "lhs";
1014 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
1015 input Boolean expand;
1016 output DAE.Exp outLhs;
1017 output DAE.Exp outRhs;
1018
1019 protected
1020 list<DAE.Exp> lhs, rhs;
1021 DAE.Exp tmpLhs, e1;
1022 Boolean b;
1023 DAE.ComponentRef cr;
1024 algorithm
1025
1026
2/2
✓ Branch 1 taken 180230 times.
✓ Branch 2 taken 181961 times.
 362191 if expHasCref(inExp1, inExp3) then
1027
1028
2/2
✓ Branch 0 taken 180167 times.
✓ Branch 1 taken 63 times.
 180230 if expand then
1029 180167 (cr, b) := Expression.expOrDerCref(inExp3);
1030
2/2
✓ Branch 0 taken 5052 times.
✓ Branch 1 taken 175115 times.
 180167 if b then
1031 5052 (lhs, rhs) := Expression.allTermsForCref(inExp1, cr, Expression.expHasDerCref);
1032 else
1033 175115 (lhs, rhs) := Expression.allTermsForCref(inExp1, cr, Expression.expHasCrefNoPreOrStart);
1034 end if;
1035 else
1036 63 (lhs, rhs) := List.split1OnTrue(Expression.terms(inExp1), expHasCref, inExp3);
1037 end if;
1038
1039 // sort
1040 // a*f(x)*b -> c*f(x)
1041 outLhs := DAE.RCONST(0.0);
1042 tmpLhs := DAE.RCONST(0.0);
1043
2/2
✓ Branch 0 taken 183460 times.
✓ Branch 1 taken 180230 times.
 363690 for e in lhs loop
1044
2/2
✓ Branch 1 taken 10645 times.
✓ Branch 2 taken 172815 times.
 183460 if Expression.isNegativeUnary(e) then
1045
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10645 times.
 10645 DAE.UNARY(exp = e1) := e;
1046 10645 tmpLhs := expAddX(e1, tmpLhs, inExp3); // special add
1047 else
1048 172815 outLhs := expAddX(e, outLhs, inExp3); // special add
1049 end if;
1050 end for;
1051 180230 outLhs := expAddX(outLhs, Expression.negate(tmpLhs), inExp3);
1052 //rhs
1053 180230 outRhs := Expression.makeSum1(rhs);
1054 180230 (outRhs,_) := ExpressionSimplify.simplify1(outRhs);
1055 180230 (outLhs,_) := ExpressionSimplify.simplify1(outLhs);
1056
1057 else
1058 outLhs := DAE.RCONST(0.0);
1059 outRhs := inExp1;
1060 end if;
1061
1062 end preprocessingSolve5;
1063
1064 protected function unifyFunCalls
1065 "
1066 e.g.
1067 semiLinear() -> if
1068 author: Vitalij Ruge
1069 "
1070 input DAE.Exp inExp1 "lhs";
1071 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
1072 output DAE.Exp oExp;
1073 output Boolean newX;
1074 algorithm
1075 35531 (oExp,_) := Expression.traverseExpTopDown(inExp1, unifyFunCallsWork, (inExp3));
1076 35531 newX := ExpressionBasics.expEqual(oExp, inExp1);
1077 end unifyFunCalls;
1078
1079 protected function unifyFunCallsWork
1080 input DAE.Exp inExp;
1081 input DAE.Exp iT;
1082 output DAE.Exp outExp;
1083 output Boolean cont;
1084 output DAE.Exp oT;
1085 algorithm
1086 (outExp,cont,oT) := match(inExp, iT)
1087 local
1088 DAE.Exp e, e1,e2,e3, X;
1089 DAE.Type tp;
1090 /*
1091 case(DAE.CALL(path = Absyn.IDENT(name = "smooth"), expLst = {_, e}),X)
1092 guard expHasCref(e, X)
1093 then (e, true, iT);
1094
1095 case(DAE.CALL(path = Absyn.IDENT(name = "noEvent"), expLst = {e}),X)
1096 guard expHasCref(e, X)
1097 then (e, true, iT);
1098 */
1099 case(DAE.CALL(path = Absyn.IDENT(name = "semiLinear"),expLst = {e1, e2, e3}),_)
1100 guard not Expression.isZero(e1)
1101 algorithm
1102 21 tp := Expression.typeof(e1);
1103 21 e := DAE.IFEXP(DAE.RELATION(e1,DAE.GREATEREQ(tp), Expression.makeConstZero(tp),-1,NONE()),Expression.expMul(e1,e2), Expression.expMul(e1,e3));
1104 then (e,true, iT);
1105
1106 // df_der(x) = (x-old(x))/dt
1107 case(DAE.CALL(path = Absyn.IDENT(name = "$_DF$DER"),expLst = {e1}),X)
1108 guard expHasCref(e1, X)
1109 algorithm
1110 ✗ tp := Expression.typeof(e1);
1111 ✗ e2 := Expression.crefExp(ComponentReferenceBasics.makeCrefIdent(BackendDAE.symSolverDT, DAE.T_REAL_DEFAULT, {}));
1112 ✗ e3 := Expression.makePureBuiltinCall("pre", {e1}, tp);
1113 ✗ e3 := Expression.expSub(e1,e3);
1114 ✗ e := Expression.expDiv(e3,e2);
1115 then (e,true, iT);
1116
1117 else (inExp, true, iT);
1118 end match;
1119
1120 end unifyFunCallsWork;
1121
1122
1123 protected function solveFunCalls
1124 "
1125 - inline modelica functions
1126 - TODO: support annotation inverse
1127 author: Vitalij Ruge
1128 "
1129 input DAE.Exp inExp1 "lhs";
1130 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
1131 input Option<AvlTreePathFunction.Tree> functions;
1132 output DAE.Exp x;
1133 output Boolean con;
1134 algorithm
1135 (x,con) := matchcontinue inExp1
1136 local DAE.Exp funX; Boolean b;
1137 case _
1138 algorithm
1139 507 (funX,_) := Expression.traverseExpTopDown(inExp1, inlineCallX, (inExp3, functions));
1140 507 b := not ExpressionBasics.expEqual(funX, inExp1);
1141 then (funX, b);
1142 else (inExp1, false);
1143 end matchcontinue;
1144 end solveFunCalls;
1145
1146 protected function removeSimpleCalls
1147 "
1148 helprer function for preprocessingSolve
1149
1150 solve e.g.
1151 exp(x) = y
1152 log(x) = y
1153 "
1154 input DAE.Exp inExp1 "lhs";
1155 input DAE.Exp inExp2 "rhs";
1156 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
1157
1158 output DAE.Exp outLhs;
1159 output DAE.Exp outRhs;
1160 output Boolean con "continue";
1161 algorithm
1162 (outLhs, outRhs, con) := match inExp1
1163 9224 case DAE.CALL() then removeSimpleCalls2(inExp1, inExp2, inExp3);
1164 5586 else (inExp1, inExp2, false);
1165 end match;
1166 end removeSimpleCalls;
1167
1168
1169 protected function removeSimpleCalls2
1170 "
1171 helprer function for preprocessingSolve
1172
1173 solve e.g.
1174 exp(x) = y
1175 log(x) = y
1176 "
1177 input DAE.Exp inExp1 "lhs";
1178 input DAE.Exp inExp2 "rhs";
1179 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
1180
1181 output DAE.Exp outLhs;
1182 output DAE.Exp outRhs;
1183 output Boolean con "continue";
1184 algorithm
1185 (outLhs, outRhs, con) := matchcontinue (inExp1, inExp2)
1186 local
1187 DAE.Exp e1, e2, e3;
1188
1189
1190 //tanh(x) =y -> x = 1/2 * ln((1+y)/(1-y))
1191 case (DAE.CALL(path = Absyn.IDENT(name = "tanh"),expLst = {e1}), _)
1192 algorithm
1193
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
 12 true := expHasCref(e1, inExp3);
1194
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
 12 false := expHasCref(inExp2, inExp3);
1195
2/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
 12 true := not(Expression.isCref(inExp2) or Expression.isConst(inExp2));
1196 12 e2 := Expression.expAdd(DAE.RCONST(1.0), inExp2);
1197 12 e3 := Expression.expSub(DAE.RCONST(1.0), inExp2);
1198 12 e2 := Expression.makeDiv(e2, e3);
1199 12 e2 := Expression.makePureBuiltinCall("log",{e2},DAE.T_REAL_DEFAULT);
1200 12 e2 := Expression.expMul(DAE.RCONST(0.5), e2);
1201 then (e1, e2, true);
1202 // sinh(x) -> ln(y+(sqrt(1+y^2))
1203 case (DAE.CALL(path = Absyn.IDENT(name = "sinh"),expLst = {e1}), _)
1204 algorithm
1205 ✗ true := expHasCref(e1, inExp3);
1206 ✗ false := expHasCref(inExp2, inExp3);
1207 ✗ true := not(Expression.isCref(inExp2) or Expression.isConst(inExp2));
1208 ✗ e2 := Expression.expPow(inExp2, DAE.RCONST(2.0));
1209 ✗ e3 := Expression.expAdd(e2,DAE.RCONST(1.0));
1210 ✗ e2 := Expression.makePureBuiltinCall("sqrt",{e3},DAE.T_REAL_DEFAULT);
1211 ✗ e3 := Expression.expAdd(inExp2, e2);
1212 ✗ e2 := Expression.makePureBuiltinCall("log",{e3},DAE.T_REAL_DEFAULT);
1213 then (e1,e2,true);
1214
1215 // log10(f(a)) = g(b) => f(a) = 10^(g(b))
1216 case (DAE.CALL(path = Absyn.IDENT(name = "log10"),expLst = {e1}), _)
1217 algorithm
1218
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 6 times.
 6 true := expHasCref(e1, inExp3);
1219
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 6 times.
 6 false := expHasCref(inExp2, inExp3);
1220 6 e2 := Expression.expPow(DAE.RCONST(10.0), inExp2);
1221 then (e1, e2, true);
1222 // log(f(a)) = g(b) => f(a) = exp(g(b))
1223 case (DAE.CALL(path = Absyn.IDENT(name = "log"),expLst = {e1}), _)
1224 algorithm
1225
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 43 times.
 43 true := expHasCref(e1, inExp3);
1226
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 43 times.
 43 false := expHasCref(inExp2, inExp3);
1227 43 e2 := Expression.makePureBuiltinCall("exp",{inExp2},DAE.T_REAL_DEFAULT);
1228 then (e1, e2, true);
1229 // exp(f(a)) = g(b) => f(a) = log(g(b))
1230 case (DAE.CALL(path = Absyn.IDENT(name = "exp"),expLst = {e1}), _)
1231 algorithm
1232
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 30 times.
 30 true := expHasCref(e1, inExp3);
1233
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 30 times.
 30 false := expHasCref(inExp2, inExp3);
1234 30 e2 := Expression.makePureBuiltinCall("log",{inExp2},DAE.T_REAL_DEFAULT);
1235 then (e1, e2, true);
1236 // sqrt(f(a)) = g(b) => f(a) = (g(b))^2
1237 case (DAE.CALL(path = Absyn.IDENT(name = "sqrt"),expLst = {e1}), _)
1238 algorithm
1239
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 26 times.
 26 true := expHasCref(e1, inExp3);
1240
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 26 times.
 26 false := expHasCref(inExp2, inExp3);
1241 e2 := DAE.RCONST(2.0);
1242 26 e2 := Expression.expPow(inExp2,e2);
1243 then (e1, e2, true);
1244 // semiLinear(0, a, b) = 0 => a = b // rule 1
1245 case (DAE.CALL(path = Absyn.IDENT(name = "semiLinear"),expLst = {DAE.RCONST(real = 0.0), e1, e2}), DAE.RCONST(real = 0.0))
1246 then (e1,e2,true);
1247 // smooth(i,f(a)) = rhs -> f(a) = rhs
1248 //case (DAE.CALL(path = Absyn.IDENT(name = "smooth"),expLst = {_, e2}),_,_)
1249 // then (e2, inExp2, true);
1250 // noEvent(f(a)) = rhs -> f(a) = rhs
1251 //case (DAE.CALL(path = Absyn.IDENT(name = "noEvent"),expLst = {e2}),_,_)
1252 // then (e2, inExp2, true);
1253
1254 else (inExp1, inExp2, false);
1255 end matchcontinue;
1256 end removeSimpleCalls2;
1257
1258 protected function inlineCallX
1259 "
1260 inline function call if depends on X where X is cref or der(cref)
1261 DAE.Exp inExp2 DAE.CREF or 'der(DAE.CREF())'
1262 author: vitalij
1263 "
1264 input DAE.Exp inExp;
1265 input tuple<DAE.Exp, Option<AvlTreePathFunction.Tree>> iT;
1266 output DAE.Exp outExp;
1267 output Boolean cont;
1268 output tuple<DAE.Exp, Option<AvlTreePathFunction.Tree>> oT;
1269 algorithm
1270 (outExp,cont,oT) := matchcontinue(inExp, iT)
1271 local
1272 DAE.Exp e, X;
1273 Option<AvlTreePathFunction.Tree> functions;
1274 Boolean b;
1275
1276 case(DAE.CALL(),(X, functions))
1277 guard expHasCref(inExp, X)
1278 algorithm
1279 //print("\nfIn: ");print(ExpressionBasics.printExpStr(inExp));
1280 372 (e,_,b) := Inline.forceInlineExp(inExp,(functions,{DAE.NORM_INLINE(),DAE.DEFAULT_INLINE()}),DAE.emptyElementSource,Ceval.cevalSimpleWithFunctionTreeReturnExp);
1281 //print("\nfOut: ");print(ExpressionBasics.printExpStr(e));
1282 372 then (e, not b, iT);
1283 else (inExp, true, iT);
1284 end matchcontinue;
1285 end inlineCallX;
1286
1287 protected function preprocessingSolveTmpVars
1288 "
1289 helper function for solveWork
1290 creat tmp vars if needed!
1291 e.g. for solve abs
1292 "
1293 input DAE.Exp inExp1;
1294 input DAE.Exp inExp2;
1295 input DAE.Exp inExp3;
1296 input Option<DAE.Exp> optCond "condition from an if expression";
1297 input Integer uniqueEqIndex "offset for tmp vars";
1298 input list<BackendDAE.Equation> ieqnForNewVars;
1299 input list<DAE.ComponentRef> inewVarsCrefs;
1300 input Integer idepth "depth of tmp var";
1301 output DAE.Exp x;
1302 output DAE.Exp y;
1303 output Boolean new_x;
1304 output list<BackendDAE.Equation> eqnForNewVars "eqn for tmp vars";
1305 output list<DAE.ComponentRef> newVarsCrefs;
1306 output Integer odepth;
1307 algorithm
1308 (x, y, new_x, eqnForNewVars, newVarsCrefs, odepth) := matchcontinue inExp1
1309 local
1310 DAE.Exp arg, e1, e_1, e2, exP, lhs, e3, e5, e6, rhs, a1,x1, a2,x2, ee1, ee2;
1311 list<DAE.Exp> z1, z2, z3, z4;
1312 DAE.Type tp;
1313 list<BackendDAE.Equation> eqnForNewVars_;
1314 list<DAE.ComponentRef> newVarsCrefs_;
1315 DAE.Operator op1;
1316 String name;
1317
1318 // try to invert a function call
1319 // f(x) = y -> x = f^(-1)(y)
1320 case DAE.CALL(path = Absyn.IDENT(name = name),expLst = {arg})
1321 guard expHasCref(arg, inExp3) and not expHasCref(inExp2, inExp3)
1322 algorithm
1323 103 (y, new_x, eqnForNewVars_, newVarsCrefs_, odepth) := preprocessingSolveFunctionCall(name, arg, inExp2, inExp3, optCond, uniqueEqIndex, idepth);
1324
1325
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 103 times.
 103 if listEmpty(eqnForNewVars_) then
1326 ✗ eqnForNewVars_ := ieqnForNewVars;
1327 else
1328 eqnForNewVars_ := match optCond
1329 local
1330 DAE.Exp cond;
1331 case SOME(cond)
1332 8 then BackendDAE.IF_EQUATION({cond}, {eqnForNewVars_}, {}, DAE.emptyElementSource, BackendDAE.EQ_ATTR_DEFAULT_UNKNOWN) :: ieqnForNewVars;
1333 95 else listAppend(eqnForNewVars_, ieqnForNewVars);
1334 end match;
1335 end if;
1336
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 103 times.
 103 then (if new_x then arg else inExp1, y, new_x, eqnForNewVars_, listAppend(newVarsCrefs_, inewVarsCrefs), odepth);
1337
1338 // x^n = y -> x = y^(1/n)
1339 case DAE.BINARY(e1, DAE.POW(tp), e2)
1340 guard expHasCref(e1, inExp3) and not expHasCref(e2, inExp3)
1341 algorithm
1342 100 tp := Expression.typeof(e1);
1343 100 exP := makeInitialGuess(tp, inExp3, e1);
1344 // exP := makeInitialGuess(tp, inExp3, inExp2);
1345 100 (exP, eqnForNewVars_, newVarsCrefs_) := makeTmpEqnAndCrefFromExp(exP, tp, "X$ABS", uniqueEqIndex, idepth, ieqnForNewVars, inewVarsCrefs, false);
1346 100 e_1 := Expression.makePureBuiltinCall("$_signNoNull", {exP}, tp); // sign
1347
1348 100 lhs := Expression.expPow(inExp2, Expression.inverseFactors(e2)); // y^(1/n)
1349 100 lhs := Expression.makePureBuiltinCall("abs", {lhs}, tp); // abs(y^(1/n))
1350 //lhs := Expression.makePureBuiltinCall("abs", {inExp2}, tp); // abs(y)
1351 //lhs := Expression.expPow(lhs, Expression.inverseFactors(e2)); // abs(y)^(1/n)
1352 100 lhs := Expression.expMul(e_1, lhs); // sign*abs(y^(1/n))
1353 100 then(e1, lhs, true, eqnForNewVars_, newVarsCrefs_, idepth + 1);
1354
1355 //QE
1356 // a*x^n + b*x^m = c
1357 // a*x^n - b*x^m = c
1358 case DAE.BINARY(ee1, op1, ee2)
1359 guard(Expression.isAddOrSub(op1))
1360 algorithm
1361 704 (z1, z2) := List.split1OnTrue(Expression.factors(ee1), expHasCref, inExp3);
1362 704 (z3, z4) := List.split1OnTrue(Expression.factors(ee2), expHasCref, inExp3);
1363
1364 704 x1 := makeProductLstSort(z1);
1365 704 a1 := makeProductLstSort(z2);
1366
1367 704 x2 := makeProductLstSort(z3);
1368
2/2
✓ Branch 1 taken 224 times.
✓ Branch 2 taken 480 times.
 704 a2 := if Expression.isAdd(op1) then makeProductLstSort(z4) else Expression.negate(makeProductLstSort(z4));
1369 704 (e2, e3) := simplifyBinaryMulCoeff(x1);
1370 704 (e5, e6) := simplifyBinaryMulCoeff(x2);
1371 704 (lhs, rhs, eqnForNewVars_, newVarsCrefs_) := solveQE(a1,e2,e3,a2,e5,e6,inExp2,inExp3,ieqnForNewVars,inewVarsCrefs,uniqueEqIndex,idepth);
1372 56 then(lhs, rhs, true, eqnForNewVars_, newVarsCrefs_, idepth + 1);
1373
1374 else (inExp1, inExp2, false, ieqnForNewVars, inewVarsCrefs, idepth);
1375 end matchcontinue;
1376 end preprocessingSolveTmpVars;
1377
1378 protected function preprocessingSolveFunctionCall
1379 input String name "of the function";
1380 input DAE.Exp arg "ument of the function";
1381 input DAE.Exp rhs;
1382 input DAE.Exp inExp3 "solve for this";
1383 input Option<DAE.Exp> optCond "condition from an if expression";
1384 input Integer uniqueEqIndex "offset for tmp vars";
1385 input Integer idepth "depth of tmp var";
1386 output DAE.Exp result;
1387 output Boolean new_x;
1388 output list<BackendDAE.Equation> newEqns "eqns for tmp vars";
1389 output list<DAE.ComponentRef> newVars "tmp vars";
1390 output Integer odepth;
1391 algorithm
1392 (result, new_x, newEqns, newVars, odepth) := match name
1393 local
1394 DAE.Exp y, exP, sgn, inv, e1, e2, k1, k2, x1, x2;
1395 DAE.Type tp;
1396 list<BackendDAE.Equation> eqns;
1397 list<DAE.ComponentRef> vars;
1398 BackendDAE.Equation ass "assertion for finding domain violations";
1399
1400 // tanh(x) -> 0.5 * ln((1 + y)/(1 - y))
1401 // exists for y in (-1, 1)
1402 // unique
1403 case "tanh" algorithm
1404 ✗ tp := Expression.typeof(rhs);
1405 ✗ (y, eqns, vars) := makeTmpEqnAndCrefFromExp(rhs, tp, "Y$TANH", uniqueEqIndex, idepth, {}, {}, false);
1406 ✗ e1 := Expression.expAdd(DAE.RCONST(1.0), y); // 1 + y
1407 ✗ e2 := Expression.expSub(DAE.RCONST(1.0), y); // 1 - y
1408 ✗ e1 := Expression.makeDiv(e1, e2); // (1 + y)/(1 - y)
1409 ✗ e1 := Expression.makePureBuiltinCall("log", {e1}, tp); // ln((1 + y)/(1 - y))
1410 ✗ inv := Expression.expMul(DAE.RCONST(0.5), e1); // 0.5 * ln((1 + y)/(1 - y))
1411 ✗ ass := makeDomainAssert(name, rhs, SOME((-1.0, false)), SOME((1.0, false))); // y in (-1, 1)
1412 ✗ then (inv, true, ass :: eqns, vars, idepth + 1);
1413
1414 // sinh(x) -> ln(y + sqrt(y^2 + 1))
1415 // exixts always
1416 // unique
1417 case "sinh" algorithm
1418 ✗ tp := Expression.typeof(rhs);
1419 ✗ (y, eqns, vars) := makeTmpEqnAndCrefFromExp(rhs, tp, "Y$SINH", uniqueEqIndex, idepth, {}, {}, false);
1420 ✗ e1 := Expression.expPow(y, DAE.RCONST(2.0)); // y^2
1421 ✗ e1 := Expression.expAdd(e1, DAE.RCONST(1.0)); // y^2 + 1
1422 ✗ e1 := Expression.makePureBuiltinCall("sqrt", {e1}, tp); // sqrt(y^2 + 1)
1423 ✗ e1 := Expression.expAdd(y, e1); // y + sqrt(y^2 + 1)
1424 ✗ e1 := Expression.makePureBuiltinCall("log", {e1}, tp); // ln(y + sqrt(y^2 + 1))
1425 ✗ then (e1, true, eqns, vars, idepth + 1);
1426
1427 // cosh(x) -> ln(y + sign*sqrt(y^2 - 1))
1428 // exists for y in [1, inf)
1429 // two values, sign is -1 or 1
1430 case "cosh" algorithm
1431 4 tp := Expression.typeof(rhs);
1432 4 (y, eqns, vars) := makeTmpEqnAndCrefFromExp(rhs, tp, "Y$COSH", uniqueEqIndex, idepth, {}, {}, false);
1433
1434 4 exP := makeInitialGuess(tp, inExp3, arg);
1435 4 (exP, eqns, vars) := makeTmpEqnAndCrefFromExp(exP, tp, "SIGN$COSH", uniqueEqIndex, idepth, eqns, vars, false);
1436 4 sgn := Expression.makePureBuiltinCall("$_signNoNull", {exP}, tp); // sign
1437
1438 4 e1 := Expression.expPow(y, DAE.RCONST(2.0)); // y^2
1439 4 e1 := Expression.expSub(e1, DAE.RCONST(1.0)); // y^2 - 1
1440 4 e1 := Expression.makePureBuiltinCall("sqrt", {e1}, tp); // sqrt(y^2 - 1)
1441 4 e1 := Expression.expMul(sgn, e1); // sign*sqrt(y^2 - 1)
1442 4 e1 := Expression.expAdd(y, e1); // y + sign*sqrt(y^2 - 1)
1443 4 e1 := Expression.makePureBuiltinCall("log", {e1}, tp); // ln(y + sign*sqrt(y^2 - 1))
1444
1445 4 ass := makeDomainAssert(name, rhs, SOME((1.0, true)), NONE()); // y in [1, inf)
1446 4 then (e1, true, ass :: eqns, vars, idepth + 1);
1447
1448 // cos(x) -> sign*acos(y) + 2*pi*k
1449 // exists for y in [-1, 1]
1450 // infinitely many values, k is integer, sign is -1 or 1
1451 case "cos" algorithm
1452 33 tp := Expression.typeof(rhs);
1453 33 (y, eqns, vars) := makeTmpEqnAndCrefFromExp(rhs, tp, "Y$COS", uniqueEqIndex, idepth, {}, {}, false);
1454
1455 33 inv := Expression.makePureBuiltinCall("acos", {y}, tp);
1456 33 (inv, eqns, vars) := makeTmpEqnAndCrefFromExp(inv, tp, "INV$COS", uniqueEqIndex, idepth, eqns, vars, false);
1457
1458 33 exP := makeInitialGuess(tp, inExp3, arg);
1459 33 (exP, eqns, vars) := makeTmpEqnAndCrefFromExp(exP, tp, "PREX$COS", uniqueEqIndex, idepth, eqns, vars, false);
1460
1461 33 k1 := helpInvCos(inv, exP, tp, true);
1462 33 k2 := helpInvCos(inv, exP, tp, false);
1463 33 (k1, eqns, vars) := makeTmpEqnAndCrefFromExp(k1, tp, "k1$COS", uniqueEqIndex, idepth, eqns, vars, false);
1464 33 (k2, eqns, vars) := makeTmpEqnAndCrefFromExp(k2, tp, "k2$COS", uniqueEqIndex, idepth, eqns, vars, false);
1465
1466 33 x1 := helpInvCos2(k1, inv, tp, true);
1467 33 x2 := helpInvCos2(k2, inv, tp, false);
1468 33 (x1, eqns, vars) := makeTmpEqnAndCrefFromExp(x1, tp, "x1$COS", uniqueEqIndex, idepth, eqns, vars, false);
1469 33 (x2, eqns, vars) := makeTmpEqnAndCrefFromExp(x2, tp, "x2$COS", uniqueEqIndex, idepth, eqns, vars, false);
1470 33 e1 := helpInvCos3(x1, x2, exP, tp);
1471
1472 33 ass := makeDomainAssert(name, rhs, SOME((-1.0, true)), SOME((1.0, true))); // y in [-1, 1]
1473 33 then (e1, true, ass :: eqns, vars, idepth + 1);
1474
1475 // sin(x) -> (-1)^k * asin(y) + k*pi
1476 // exists for y in [-1, 1]
1477 // infinitely many values, k is integer
1478 case "sin" algorithm
1479 27 tp := Expression.typeof(rhs);
1480 27 (y, eqns, vars) := makeTmpEqnAndCrefFromExp(rhs, tp, "Y$SIN", uniqueEqIndex, idepth, {}, {}, false);
1481
1482 27 inv := Expression.makePureBuiltinCall("asin", {y}, tp);
1483 27 (inv, eqns, vars) := makeTmpEqnAndCrefFromExp(inv, tp, "INV$SIN", uniqueEqIndex, idepth, eqns, vars, false);
1484
1485 27 exP := makeInitialGuess(tp,inExp3,arg);
1486 27 (exP, eqns, vars) := makeTmpEqnAndCrefFromExp(exP, tp, "PREX$SIN", uniqueEqIndex, idepth, eqns, vars, false);
1487
1488 27 k1 := helpInvSin(inv, arg, tp, true);
1489 27 k2 := helpInvSin(inv, arg, tp, false);
1490 27 (k1, eqns, vars) := makeTmpEqnAndCrefFromExp(k1, tp, "k1$SIN", uniqueEqIndex, idepth, eqns, vars, false);
1491 27 (k2, eqns, vars) := makeTmpEqnAndCrefFromExp(k2, tp, "k2$SIN", uniqueEqIndex, idepth, eqns, vars, false);
1492
1493 27 x1 := helpInvSin2(k1, inv, tp, true);
1494 27 x2 := helpInvSin2(k2, inv, tp, false);
1495 27 (x1, eqns, vars) := makeTmpEqnAndCrefFromExp(x1, tp, "x1$SIN", uniqueEqIndex, idepth, eqns, vars, false);
1496 27 (x2, eqns, vars) := makeTmpEqnAndCrefFromExp(x2, tp, "x2$SIN", uniqueEqIndex, idepth, eqns, vars, false);
1497
1498 27 e1 := helpInvCos3(x1, x2, exP, tp);
1499
1500 27 ass := makeDomainAssert(name, rhs, SOME((-1.0, true)), SOME((1.0, true))); // y in [-1, 1]
1501 27 then (e1, true, ass :: eqns, vars, idepth + 1);
1502
1503 // tan(x) -> atan(y) + k*pi
1504 // exists always
1505 // infinitely many values, k is integer
1506 case "tan" algorithm
1507 5 tp := Expression.typeof(rhs);
1508 5 (y, eqns, vars) := makeTmpEqnAndCrefFromExp(rhs, tp, "Y$TAN", uniqueEqIndex, idepth, {}, {}, false);
1509
1510 5 inv := Expression.makePureBuiltinCall("atan", {y}, tp);
1511 5 (inv, eqns, vars) := makeTmpEqnAndCrefFromExp(inv, tp, "INV$TAN", uniqueEqIndex, idepth, eqns, vars, false);
1512
1513 5 exP := makeInitialGuess(tp, inExp3, arg);
1514 5 (exP, eqns, vars) := makeTmpEqnAndCrefFromExp(exP, tp, "PREX$TAN", uniqueEqIndex, idepth, eqns, vars, false);
1515
1516 5 k1 := Expression.expSub(exP, inv); // pre(x) - atan(y)
1517 5 k1 := Expression.makeDiv(k1, DAE.PI); // (pre(x) - atan(y))/pi
1518 5 k1 := Expression.makePureBuiltinCall("$_round", {k1}, tp); // k = round((pre(x) - atan(y))/pi)
1519 5 e1 := Expression.expMul(k1, DAE.PI); // k*pi
1520 5 e1 := Expression.expAdd(inv, e1); // atan(y) + pi*k
1521 5 then (e1, true, eqns, vars, idepth + 1);
1522
1523 // abs(x) -> sign*y
1524 // exists for y in [0, inf)
1525 // two values, sign is -1 or 1
1526 case "abs" algorithm
1527 34 tp := Expression.typeof(arg);
1528 34 exP := makeInitialGuess(tp, inExp3, arg);
1529 34 (exP, eqns, vars) := makeTmpEqnAndCrefFromExp(exP, tp, "SIGN$ABS", uniqueEqIndex, idepth, {}, {}, false);
1530 34 sgn := Expression.makePureBuiltinCall("$_signNoNull", {exP}, tp); // sign
1531 34 e1 := Expression.expMul(sgn, rhs); // sign*y
1532 34 ass := makeDomainAssert(name, rhs, SOME((0.0, true)), NONE()); // y in [0, inf)
1533 34 then (e1, true, ass ::eqns, vars, idepth + 1);
1534
1535 // sqrt(x) -> y^2
1536 // exists for y in [0, inf)
1537 // unique
1538 case "sqrt" algorithm
1539 ✗ inv := Expression.expPow(rhs, DAE.RCONST(2.0)); // y^2
1540 ✗ ass := makeDomainAssert(name, rhs, SOME((0.0, true)), NONE()); // y in [0, inf)
1541 ✗ then (inv, true, {ass}, {}, idepth + 1);
1542
1543 // asin(x) -> sin(y)
1544 // exists for y in [-pi/2, pi/2]
1545 // unique
1546 case "asin" algorithm
1547 ✗ tp := Expression.typeof(rhs);
1548 ✗ inv := Expression.makePureBuiltinCall("sin", {rhs}, tp); // sin(y)
1549 ✗ ass := makeDomainAssert(name, rhs, SOME((-0.5*Expression.toReal(DAE.PI), true)), SOME((0.5*Expression.toReal(DAE.PI), true))); // y in [-pi/2, pi/2]
1550 ✗ then (inv, true, {ass}, {}, idepth + 1);
1551
1552 // acos(x) -> cos(y)
1553 // exists for y in [0, pi]
1554 // unique
1555 case "acos" algorithm
1556 ✗ tp := Expression.typeof(rhs);
1557 ✗ inv := Expression.makePureBuiltinCall("cos", {rhs}, tp); // cos(y)
1558 ✗ ass := makeDomainAssert(name, rhs, SOME((0.0, true)), SOME((Expression.toReal(DAE.PI), true))); // y in [0, pi]
1559 ✗ then (inv, true, {ass}, {}, idepth + 1);
1560
1561 // atan(x) -> tan(y)
1562 // exists for y in [-pi/2, pi/2]
1563 // unique
1564 case "atan" algorithm
1565 ✗ tp := Expression.typeof(rhs);
1566 ✗ inv := Expression.makePureBuiltinCall("tan", {rhs}, tp); // tan(y)
1567 ✗ ass := makeDomainAssert(name, rhs, SOME((-0.5*Expression.toReal(DAE.PI), true)), SOME((0.5*Expression.toReal(DAE.PI), true))); // y in [-pi/2, pi/2]
1568 ✗ then (inv, true, {ass}, {}, idepth + 1);
1569
1570 // exp(x) -> log(y)
1571 // exists for y in (0, inf)
1572 // unique
1573 case "exp" algorithm
1574 ✗ tp := Expression.typeof(rhs);
1575 ✗ inv := Expression.makePureBuiltinCall("log", {rhs}, tp); // log(y)
1576 ✗ ass := makeDomainAssert(name, rhs, SOME((0.0, false)), NONE()); // y in (0, inf)
1577 ✗ then (inv, true, {ass}, {}, idepth + 1);
1578
1579 // log(x) -> exp(y)
1580 // exists always
1581 // unique
1582 case "log" algorithm
1583 ✗ tp := Expression.typeof(rhs);
1584 ✗ inv := Expression.makePureBuiltinCall("exp", {rhs}, tp); // exp(y)
1585 ✗ then (inv, true, {}, {}, idepth + 1);
1586
1587 // log10(x) -> 10^y
1588 // exists always
1589 // unique
1590 case "log10" algorithm
1591 ✗ inv := Expression.expPow(DAE.RCONST(10.0), rhs); // 10^y
1592 ✗ then (inv, true, {}, {}, idepth + 1);
1593
1594 // sign(x) is not invertible
1595 case "sign" then (rhs, false, {}, {}, idepth);
1596
1597 // $_DF$DER(x) = y -> (x - pre(x))/dt = y -> x = y*dt + pre(x)
1598 case "$_DF$DER" algorithm
1599 ✗ e1 := Expression.crefExp(ComponentReferenceBasics.makeCrefIdent(BackendDAE.symSolverDT, DAE.T_REAL_DEFAULT, {})); // dt
1600 ✗ exP := Expression.makePureBuiltinCall("pre", {arg}, Expression.typeof(arg)); // pre(x)
1601 ✗ e1 := Expression.expAdd(Expression.expMul(rhs, e1), exP); // y*dt + pre(x)
1602 ✗ then(e1, true, {}, {}, idepth + 1);
1603
1604 // don't know inverse of this function
1605 else (rhs, false, {}, {}, idepth);
1606 end match;
1607 end preprocessingSolveFunctionCall;
1608
1609 protected function simplifyBinaryMulCoeff
1610 "generalization of ExpressionSimplify.simplifyBinaryMulCoeff2"
1611 input DAE.Exp inExp;
1612 output DAE.Exp exp1;
1613 output DAE.Exp exp2;
1614 algorithm
1615 (exp1, exp2) := match inExp
1616 local
1617 DAE.Exp e,e1,e2;
1618 DAE.Exp coeff;
1619
1620 case e as DAE.CREF()
1621 then ((e, DAE.RCONST(1.0)));
1622
1623 case DAE.BINARY(exp1 = e1,operator = DAE.POW(),exp2 = DAE.UNARY(operator = DAE.UMINUS(), exp = coeff))
1624 ✗ then
1625 ((e1, Expression.negate(coeff)));
1626
1627 case DAE.BINARY(exp1 = e1,operator = DAE.POW(),exp2 = coeff)
1628 then ((e1,coeff));
1629
1630 case DAE.BINARY(exp1 = e1,operator = DAE.MUL(),exp2 = e2)
1631 guard(ExpressionBasics.expEqual(e1, e2))
1632 then
1633 ((e1, DAE.RCONST(2.0)));
1634
1635 case DAE.BINARY(e1, DAE.DIV(), e2)
1636 guard(Expression.isOne(e1))
1637 then(e2, DAE.RCONST(-1.0));
1638
1639 case DAE.CALL(path=Absyn.IDENT("sqrt"),expLst={e})
1640 then ((e,DAE.RCONST(0.5)));
1641
1642 else ((inExp,DAE.RCONST(1.0)));
1643
1644 end match;
1645 end simplifyBinaryMulCoeff;
1646
1647 protected function solveQE
1648 "
1649 solve Quadratic equation with respect to inExp3
1650 IN: a,x,n,b,y,m
1651 where solve a*x^n + b*y^m = inExp2 with 2*m = n or 2*n = m and y = x
1652
1653 author: Vitalij Ruge
1654 "
1655 input DAE.Exp e1,e2,e3,e4,e5,e6;
1656 input DAE.Exp inExp2;
1657 input DAE.Exp inExp3;
1658
1659 input list<BackendDAE.Equation> ieqnForNewVars "eqn for tmp vars";
1660 input list<DAE.ComponentRef> inewVarsCrefs "cref for tmp vars";
1661 input Integer uniqueEqIndex, idepth "need for tmp vars";
1662
1663 output DAE.Exp rhs;
1664 output DAE.Exp lhs;
1665 output list<BackendDAE.Equation> eqnForNewVars;
1666 output list<DAE.ComponentRef> newVarsCrefs;
1667
1668 protected
1669 DAE.Exp e7, con, invExp, x1, x2, x, exP;
1670 DAE.Exp a,b,c, n, sgnb, b2, ac, sExp1, sExp2;
1671 DAE.Type tp;
1672 Boolean b1, b3;
1673 algorithm
1674
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 704 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 704 false := Expression.isZero(e1) and Expression.isZero(e2);
1675
2/2
✓ Branch 1 taken 644 times.
✓ Branch 2 taken 60 times.
 704 true := ExpressionBasics.expEqual(e2,e5);
1676 60 b1 := ExpressionBasics.expEqual(e3, Expression.expMul(DAE.RCONST(2.0),e6));
1677 60 b3 := ExpressionBasics.expEqual(e6, Expression.expMul(DAE.RCONST(2.0),e3));
1678
1679
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 56 times.
 60 true := b1 or b3;
1680
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 56 times.
 56 false := expHasCref(e1, inExp3);
1681
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 56 times.
 56 true := expHasCref(e2, inExp3);
1682
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 56 times.
 56 false := expHasCref(e3, inExp3);
1683
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 56 times.
 56 false := expHasCref(e4, inExp3);
1684
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 56 times.
 56 true := expHasCref(e5, inExp3);
1685
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 56 times.
 56 false := expHasCref(e6, inExp3);
1686
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 56 times.
 56 false := expHasCref(inExp2, inExp3);
1687
1688
1689
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 6 times.
 56 a := if b1 then e1 else e4;
1690
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 6 times.
 56 b := if b1 then e4 else e1;
1691 56 c := Expression.negate(inExp2);
1692
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 6 times.
 56 n := if b1 then e6 else e3;
1693
1694 56 tp := Expression.typeof(a);
1695 56 (a, eqnForNewVars, newVarsCrefs) := makeTmpEqnAndCrefFromExp(a, tp, "a$QE", uniqueEqIndex, idepth, ieqnForNewVars, inewVarsCrefs, false);
1696 56 con := DAE.RELATION(a,DAE.EQUAL(tp),DAE.RCONST(0.0),-1,NONE());
1697
1698 56 tp := Expression.typeof(b);
1699 56 (b, eqnForNewVars, newVarsCrefs) := makeTmpEqnAndCrefFromExp(b, tp, "b$QE", uniqueEqIndex, idepth, eqnForNewVars, newVarsCrefs, false);
1700 56 sgnb := Expression.makePureBuiltinCall("$_signNoNull",{b},tp);
1701 56 b2 := Expression.expPow(b, DAE.RCONST(2.0));
1702 56 (b2, eqnForNewVars, newVarsCrefs) := makeTmpEqnAndCrefFromExp(b2, tp, "bPow2$QE", uniqueEqIndex, idepth, eqnForNewVars, newVarsCrefs, false);
1703
1704 56 tp := Expression.typeof(c);
1705 56 (c, eqnForNewVars, newVarsCrefs) := makeTmpEqnAndCrefFromExp(c, tp, "c$QE", uniqueEqIndex, idepth, eqnForNewVars, newVarsCrefs, false);
1706 56 ac := Expression.expMul(a,c);
1707 56 ac := Expression.expMul(DAE.RCONST(4.0),ac);
1708
1709 56 sExp1 := Expression.expSub(b2,ac);
1710 56 sExp2 := Expression.makePureBuiltinCall("sqrt",{sExp1},tp);
1711 56 sExp2 := Expression.expMul(sgnb, sExp2);
1712
1713
1714 56 a := DAE.IFEXP(con, Expression.makeConstOne(tp), a);
1715 56 (a, eqnForNewVars, newVarsCrefs) := makeTmpEqnAndCrefFromExp(a, tp, "a1$QE", uniqueEqIndex, idepth, eqnForNewVars, newVarsCrefs, false);
1716
1717 56 x1 := Expression.expAdd(b, sExp2);
1718 56 x1 := Expression.makeDiv(x1, a);
1719 56 x1 := Expression.expMul(DAE.RCONST(-0.5), x1);
1720 56 tp := Expression.typeof(x1);
1721 56 x1 := DAE.IFEXP(con, Expression.makeConstOne(tp), x1);
1722 56 (x1, eqnForNewVars, newVarsCrefs) := makeTmpEqnAndCrefFromExp(x1, tp, "x1$QE", uniqueEqIndex, idepth, eqnForNewVars, newVarsCrefs, false);
1723
1724 //Vieta
1725 56 x2 := Expression.expMul(a,x1);
1726 56 x2 := Expression.makeDiv(c,x2);
1727 56 x2 := DAE.IFEXP(con, Expression.makeConstOne(tp), x2);
1728 56 x2 := DAE.IFEXP(DAE.RELATION(x1,DAE.EQUAL(tp),DAE.RCONST(0.0),-1,NONE()), DAE.RCONST(0.0), x2);
1729 56 (x2, eqnForNewVars, newVarsCrefs) := makeTmpEqnAndCrefFromExp(x2, tp, "x2$QE", uniqueEqIndex, idepth, eqnForNewVars, newVarsCrefs, false);
1730
1731 56 tp := Expression.typeof(e2);
1732 56 exP := makeInitialGuess(tp,inExp3,e2);
1733 56 (exP, eqnForNewVars, newVarsCrefs) := makeTmpEqnAndCrefFromExp(exP, tp, "prex$QE", uniqueEqIndex, idepth, eqnForNewVars, newVarsCrefs, false);
1734
1735 56 x := helpInvCos3(x1,x2,exP,tp);
1736 56 (x, eqnForNewVars, newVarsCrefs) := makeTmpEqnAndCrefFromExp(x, tp, "x$QE", uniqueEqIndex, idepth, eqnForNewVars, newVarsCrefs, false);
1737
1738 // a = 0
1739 56 e7 := Expression.makeDiv(inExp2,b);
1740 56 invExp := Expression.inverseFactors(n);
1741 56 (invExp, _) := ExpressionSimplify.simplify1(invExp);
1742 56 e7 := Expression.expPow(e7, invExp);
1743
1744 // if a==0
1745 56 rhs := DAE.IFEXP(con, e7 , x);
1746 // lhs
1747
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 50 times.
 56 lhs := if b1 then Expression.expPow(e2, e6) else Expression.expPow(e2, e3);
1748
1749 end solveQE;
1750
1751 protected function solveIfExp
1752 "
1753 solve:
1754 if(f(y), f(x), g(x) ) = h(y) w.r.t. x
1755 "
1756 input DAE.Exp inExp1;
1757 input DAE.Exp inExp2;
1758 input DAE.Exp inExp3;
1759 input Option<DAE.Exp> inCond;
1760 input Option<AvlTreePathFunction.Tree> functions;
1761 input Option<Integer> uniqueEqIndex "offset for tmp vars";
1762 input Integer idepth;
1763 input Boolean doInline;
1764 input Boolean isContinuousIntegration;
1765 output DAE.Exp outExp;
1766 output list<DAE.Statement> outAsserts;
1767 output list<BackendDAE.Equation> eqnForNewVars "eqn for tmp vars";
1768 output list<DAE.ComponentRef> newVarsCrefs;
1769 output Integer odepth;
1770 algorithm
1771 (outExp, outAsserts, eqnForNewVars, newVarsCrefs, odepth) := match inExp1
1772 local
1773 DAE.Exp eCond, eThen, eElse, res, lhs, rhs, cond1, cond2;
1774 list<DAE.Statement> asserts, asserts1;
1775 list<BackendDAE.Equation> eqns, eqns1;
1776 list<DAE.ComponentRef> var, var1;
1777 Integer depth;
1778
1779 // f(a) if(g(b)) then f1(a) else f2(a) =>
1780 // a1 = solve f(a),f1(a) for a
1781 // a2 = solve f(a),f2(a) for a
1782 // => a = if g(b) then a1 else a2
1783 case DAE.IFEXP(eCond, eThen, eElse)
1784 guard
1785 isContinuousIntegration or not expHasCref(eCond, inExp3)
1786 algorithm
1787
1788 // nested if expressions need to combine their conditions
1789 (cond1, cond2) := match inCond
1790 local DAE.Exp theCond;
1791 case SOME(theCond)
1792 90 then (DAE.LBINARY(theCond, DAE.AND(Expression.typeof(eCond)), eCond), DAE.LBINARY(theCond, DAE.AND(Expression.typeof(eCond)), Expression.negate(eCond)));
1793 255 else (eCond, Expression.negate(eCond));
1794 end match;
1795
1796 345 (lhs, asserts1, eqns, var, depth) := solveWork(eThen, inExp2, inExp3, SOME(cond1), functions, uniqueEqIndex, idepth, doInline, isContinuousIntegration);
1797 299 (rhs, _, eqns1, var1, depth) := solveWork(eElse, inExp2, inExp3, SOME(cond2), functions, uniqueEqIndex, depth, doInline, isContinuousIntegration);
1798
1799 247 res := DAE.IFEXP(eCond, lhs, rhs);
1800 247 asserts := listAppend(asserts1, asserts1);
1801
1/2
✓ Branch 2 taken 247 times.
✗ Branch 3 not taken.
 247 then
1802 (res, asserts, listAppend(eqns1, eqns), listAppend(var1, var), depth);
1803 else fail();
1804 end match;
1805 end solveIfExp;
1806
1807 protected function solveLinearSystem
1808 "
1809 solve linear system with newton step
1810
1811 ToDo:
1812 fixed is for ./simulation/modelica/equations/deriveToLog.mos
1813 "
1814 input DAE.Exp inExp1;
1815 input DAE.Exp inExp2;
1816 input DAE.Exp inExp3;
1817 input Option<AvlTreePathFunction.Tree> functions;
1818 input Integer idepth;
1819 output DAE.Exp outExp;
1820 output list<DAE.Statement> outAsserts;
1821 output list<BackendDAE.Equation> eqnForNewVars = {} "eqn for tmp vars";
1822 output list<DAE.ComponentRef> newVarsCrefs = {};
1823 output Integer odepth = idepth;
1824
1825
1826 algorithm
1827 (outExp,outAsserts) := match inExp3
1828 local
1829 DAE.Exp dere,e,z;
1830 DAE.ComponentRef cr;
1831 DAE.Exp rhs;
1832 DAE.Type tp;
1833 Integer i;
1834
1835 // cr = (e1-e2)/(der(e1-e2,cr))
1836 case DAE.CREF(componentRef = cr)
1837 algorithm
1838
2/2
✓ Branch 1 taken 43 times.
✓ Branch 2 taken 4561 times.
 4604 false := hasOnlyFactors(inExp1,inExp2);
1839 4561 e := Expression.expSub(inExp1,inExp2);
1840 4097 (e,_) := ExpressionSimplify.simplify1(e);
1841 //print("\ne: ");print(ExpressionBasics.printExpStr(e));
1842 4097 dere := Differentiate.differentiateExpSolve(e, cr, functions);
1843 //print("\nder(e): ");print(ExpressionBasics.printExpStr(dere));
1844 2240 (dere,_) := ExpressionSimplify.simplify(dere);
1845
2/2
✓ Branch 1 taken 1445 times.
✓ Branch 2 taken 795 times.
 2240 false := Expression.isZero(dere);
1846
2/2
✓ Branch 1 taken 689 times.
✓ Branch 2 taken 106 times.
 795 false := Expression.expHasCrefNoPreOrStart(dere, cr);
1847 106 tp := Expression.typeof(inExp3);
1848 106 (z,_) := Expression.makeZeroExpression(Expression.arrayDimension(tp));
1849 106 (e,i) := Expression.replaceExp(e, inExp3, z);
1850 // replace at least once, otherwise it's wrong
1851
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 101 times.
 106 if i < 1 then
1852 5 fail();
1853 end if;
1854 101 (e,_) := ExpressionSimplify.simplify(e);
1855 101 rhs := Expression.negate(Expression.makeDiv(e,dere));
1856 then
1857 (rhs,{});
1858
1859 else fail();
1860 end match;
1861
1862 end solveLinearSystem;
1863
1864 protected function hasOnlyFactors "help function to solve2, returns true if equation e1 == e2, has either e1 == 0 or e2 == 0 and the expression only contains
1865 factors, e.g. a*b*c = 0. In this case we can not solve the equation"
1866 input DAE.Exp e1;
1867 input DAE.Exp e2;
1868 output Boolean res;
1869 algorithm
1870 res := matchcontinue e2
1871
1872 // try normal
1873 case _
1874 algorithm
1875
2/2
✓ Branch 1 taken 3112 times.
✓ Branch 2 taken 1492 times.
 4604 true := Expression.isZero(e1);
1876 // More than two factors
1877
4/4
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1490 times.
✓ Branch 3 taken 1451 times.
✓ Branch 4 taken 39 times.
 1492 _::_::_ := Expression.factors(e2);
1878 //.. and more than two crefs
1879
3/4
✗ Branch 1 not taken.
✓ Branch 2 taken 39 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 36 times.
 39 _::_::_ := Expression.extractCrefsFromExp(e2);
1880 then
1881 true;
1882
1883 // swapped args
1884 case _
1885 algorithm
1886
2/2
✓ Branch 1 taken 2786 times.
✓ Branch 2 taken 1782 times.
 4568 true := Expression.isZero(e2);
1887
3/4
✗ Branch 1 not taken.
✓ Branch 2 taken 1782 times.
✓ Branch 3 taken 1757 times.
✓ Branch 4 taken 25 times.
 1782 _::_::_ := Expression.factors(e1);
1888
3/4
✗ Branch 1 not taken.
✓ Branch 2 taken 25 times.
✓ Branch 3 taken 18 times.
✓ Branch 4 taken 7 times.
 25 _::_::_ := Expression.extractCrefsFromExp(e1);
1889 then
1890 true;
1891
1892 else false;
1893
1894 end matchcontinue;
1895 end hasOnlyFactors;
1896
1897
1898 protected function expHasCref
1899 "
1900 helper function for solve.
1901 case distinction for
1902 DAE.CREF or 'der(DAE.CREF())'
1903 Expression.expHasCrefNoPreOrStart
1904 or
1905 Expression.expHasDerCref
1906 "
1907 input DAE.Exp inExp1;
1908 input DAE.Exp inExp3 "DAE.CREF or 'der(DAE.CREF())'";
1909 output Boolean res;
1910
1911 algorithm
1912 res := match inExp3
1913 local DAE.ComponentRef cr;
1914
1915 1234604 case DAE.CREF(componentRef = cr) then Expression.expHasCrefNoPreOrStart(inExp1, cr);
1916 38546 case DAE.CALL(path = Absyn.IDENT(name = "der"),expLst = {DAE.CREF(componentRef = cr)}) then Expression.expHasDerCref(inExp1, cr);
1917 else
1918 algorithm
1919 ✗ if Flags.isSet(Flags.FAILTRACE) then
1920 ✗ print("\n-ExpressionSolve.solve failed:");
1921 ✗ print(" with respect to: ");print(ExpressionBasics.printExpStr(inExp3));
1922 ✗ print(" not support!");
1923 ✗ print("\n");
1924 end if;
1925 ✗ then fail();
1926 end match;
1927
1928 end expHasCref;
1929
1930 protected function makeProductLstSort
1931 "Takes a list of expressions an makes a product
1932 expression multiplying all elements in the list.
1933
1934 - a*if(b,c,d) -> if(b,a*c,a*d)
1935
1936 "
1937 input list<DAE.Exp> inExpLst;
1938 output DAE.Exp outExp;
1939 protected
1940 DAE.Type tp;
1941 list<DAE.Exp> expLstDiv, expLst, expLst2;
1942 DAE.Exp e, e1, e2;
1943 DAE.Operator op;
1944 algorithm
1945
2/2
✓ Branch 0 taken 692988 times.
✓ Branch 1 taken 808580 times.
 1501568 if listEmpty(inExpLst) then
1946 outExp := DAE.RCONST(1.0);
1947 692988 return;
1948 end if;
1949
1950 808580 tp := Expression.typeof(listHead(inExpLst));
1951
1952 808580 (expLstDiv, expLst) := List.splitOnTrue(inExpLst, Expression.isDivBinary);
1953 808580 outExp := makeProductLstSort2(expLst, tp);
1954
2/2
✓ Branch 0 taken 800774 times.
✓ Branch 1 taken 7806 times.
 808580 if not listEmpty(expLstDiv) then
1955 expLst2 := {};
1956 7806 expLst := {};
1957
1958
2/2
✓ Branch 0 taken 9443 times.
✓ Branch 1 taken 7806 times.
 17249 for elem in expLstDiv loop
1959
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9443 times.
 9443 DAE.BINARY(e1,op,e2) := elem;
1960 9443 expLst := e1::expLst;
1961 expLst2 := e2::expLst2;
1962 end for;
1963
1964
1/2
✓ Branch 0 taken 7806 times.
✗ Branch 1 not taken.
 7806 if not listEmpty(expLst2) then
1965 7806 e := makeProductLstSort(expLst2);
1966
1/2
✓ Branch 1 taken 7806 times.
✗ Branch 2 not taken.
 7806 if not Expression.isOne(e) then
1967 7806 outExp := Expression.makeDiv(outExp,e);
1968 end if;
1969 end if;
1970
1971
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7806 times.
 7806 if not listEmpty(expLst) then
1972 7806 e := makeProductLstSort(expLst);
1973 7806 outExp := Expression.expMul(outExp,e);
1974 end if;
1975
1976 end if;
1977
1978 end makeProductLstSort;
1979
1980
1981 protected function makeProductLstSort2
1982 input list<DAE.Exp> inExpLst;
1983 input DAE.Type tp;
1984 output DAE.Exp outExp = Expression.makeConstOne(tp);
1985 protected
1986 list<DAE.Exp> rest;
1987 algorithm
1988 808580 rest := ExpressionSimplify.simplifyList(inExpLst);
1989
2/2
✓ Branch 0 taken 837421 times.
✓ Branch 1 taken 808580 times.
 1646001 for elem in rest loop
1990
2/2
✓ Branch 1 taken 827962 times.
✓ Branch 2 taken 9459 times.
 837421 if not Expression.isOne(elem) then
1991 outExp := match elem
1992 local DAE.Exp e1,e2,e3;
1993 case DAE.IFEXP(e1,e2,e3)
1994 1026 then DAE.IFEXP(e1, Expression.expMul(outExp,e2), Expression.expMul(outExp,e3));
1995 826936 else Expression.expMul(outExp, elem);
1996 end match;
1997 end if;
1998 end for;
1999
2000 end makeProductLstSort2;
2001
2002 protected function makeTmpEqnAndCrefFromExp
2003 input DAE.Exp iExp;
2004 input DAE.Type tp;
2005 input String name;
2006 input Integer index1, index2;
2007 input list<BackendDAE.Equation> ieqnForNewVars;
2008 input list<DAE.ComponentRef> inewVarsCrefs;
2009 input Boolean need;
2010 output DAE.Exp oExp;
2011 output list<BackendDAE.Equation> oeqnForNewVars;
2012 output list<DAE.ComponentRef> onewVarsCrefs;
2013 protected
2014 DAE.ComponentRef cr;
2015 BackendDAE.Equation eqn;
2016 algorithm
2017 1081 (oExp,_) := ExpressionSimplify.simplify1(iExp);
2018
5/6
✓ Branch 0 taken 1081 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 971 times.
✓ Branch 4 taken 110 times.
✓ Branch 6 taken 849 times.
✓ Branch 7 taken 122 times.
 1081 if need or not (Expression.isCref(oExp) or Expression.isConst(oExp)) then
2019 849 cr := ComponentReferenceBasics.makeCrefIdent("$TMP$VAR$" + intString(index1) + "$" + intString(index2) + name, tp , {});
2020 849 eqn := BackendDAE.SOLVED_EQUATION(cr, oExp, DAE.emptyElementSource, BackendDAE.EQ_ATTR_DEFAULT_UNKNOWN);
2021 849 oExp := Expression.crefExp(cr);
2022 oeqnForNewVars := eqn :: ieqnForNewVars;
2023 849 onewVarsCrefs := cr :: inewVarsCrefs;
2024 else
2025 oeqnForNewVars := ieqnForNewVars;
2026 onewVarsCrefs := inewVarsCrefs;
2027 end if;
2028 end makeTmpEqnAndCrefFromExp;
2029
2030 protected function makeDomainAssert
2031 input String name "of the function";
2032 input DAE.Exp rhs "solution of the function";
2033 input Option<tuple<Real, Boolean>> lowerBound "(value, including?)";
2034 input Option<tuple<Real, Boolean>> upperBound "(value, including?)";
2035 output BackendDAE.Equation assEq;
2036 protected
2037 String msg;
2038 DAE.Exp cond;
2039 DAE.Algorithm algo;
2040 DAE.Type tp = Expression.typeof(rhs);
2041 algorithm
2042 (msg, cond) := match (lowerBound, upperBound)
2043 local
2044 Real lower, upper;
2045 String str;
2046 DAE.Exp l, u;
2047
2048 // range [l, u]
2049 case (SOME((lower, true)), SOME((upper, true))) algorithm
2050 60 str := "Model error: Result of " + name + " outside the range "
2051 + realString(lower) + " <= " + ExpressionBasics.printExpStr(rhs)
2052 + " <= " + realString(upper) + ". Unable to invert.";
2053 60 l := DAE.RELATION(DAE.RCONST(lower), DAE.LESSEQ(tp), rhs, -1, NONE());
2054 60 u:= DAE.RELATION(rhs, DAE.LESSEQ(tp), DAE.RCONST(upper), -1, NONE());
2055 60 then (str, DAE.LBINARY(l, DAE.AND(tp), u));
2056
2057 // range [l, u)
2058 case (SOME((lower, true)), SOME((upper, false))) algorithm
2059 ✗ str := "Model error: Result of " + name + " outside the range "
2060 + realString(lower) + " <= " + ExpressionBasics.printExpStr(rhs)
2061 + " < " + realString(upper) + ". Unable to invert.";
2062 ✗ l:= DAE.RELATION(DAE.RCONST(lower), DAE.LESSEQ(tp), rhs, -1, NONE());
2063 ✗ u:= DAE.RELATION(rhs, DAE.LESS(tp), DAE.RCONST(upper), -1, NONE());
2064 ✗ then (str, DAE.LBINARY(l, DAE.AND(tp), u));
2065
2066 // range (l, u]
2067 case (SOME((lower, false)), SOME((upper, true))) algorithm
2068 ✗ str := "Model error: Result of " + name + " outside the range "
2069 + realString(lower) + " < " + ExpressionBasics.printExpStr(rhs)
2070 + " <= " + realString(upper) + ". Unable to invert.";
2071 ✗ l:= DAE.RELATION(DAE.RCONST(lower), DAE.LESS(tp), rhs, -1, NONE());
2072 ✗ u:= DAE.RELATION(rhs, DAE.LESSEQ(tp), DAE.RCONST(upper), -1, NONE());
2073 ✗ then (str, DAE.LBINARY(l, DAE.AND(tp), u));
2074
2075 // range (l, u)
2076 case (SOME((lower, false)), SOME((upper, false))) algorithm
2077 ✗ str := "Model error: Result of " + name + " outside the range "
2078 + realString(lower) + " < " + ExpressionBasics.printExpStr(rhs)
2079 + " < " + realString(upper) + ". Unable to invert.";
2080 ✗ l:= DAE.RELATION(DAE.RCONST(lower), DAE.LESS(tp), rhs, -1, NONE());
2081 ✗ u:= DAE.RELATION(rhs, DAE.LESS(tp), DAE.RCONST(upper), -1, NONE());
2082 ✗ then (str, DAE.LBINARY(l, DAE.AND(tp), u));
2083
2084 // range [l, inf)
2085 case (SOME((lower, true)), NONE()) algorithm
2086 38 str := "Model error: Result of " + name + " should be "
2087 + ExpressionBasics.printExpStr(rhs) + " >= " + realString(lower)
2088 + ". Unable to invert.";
2089 38 l:= DAE.RELATION(DAE.RCONST(lower), DAE.LESSEQ(tp), rhs, -1, NONE());
2090 then (str, l);
2091
2092 // range (l, inf)
2093 case (SOME((lower, true)), NONE()) algorithm
2094 ✗ str := "Model error: Result of " + name + " should be "
2095 + ExpressionBasics.printExpStr(rhs) + " > " + realString(lower)
2096 + ". Unable to invert.";
2097 ✗ l:= DAE.RELATION(DAE.RCONST(lower), DAE.LESS(tp), rhs, -1, NONE());
2098 then (str, l);
2099
2100 // range (-inf, u]
2101 case (NONE(), SOME((upper, true))) algorithm
2102 ✗ str := "Model error: Result of " + name + " should be "
2103 + ExpressionBasics.printExpStr(rhs) + " <= " + realString(upper)
2104 + ". Unable to invert.";
2105 ✗ u:= DAE.RELATION(rhs, DAE.LESSEQ(tp), DAE.RCONST(upper), -1, NONE());
2106 then (str, u);
2107
2108 // range (-inf, u)
2109 case (NONE(), SOME((upper, false))) algorithm
2110 ✗ str := "Model error: Result of " + name + " should be "
2111 + ExpressionBasics.printExpStr(rhs) + " < " + realString(upper)
2112 + ". Unable to invert.";
2113 ✗ u:= DAE.RELATION(rhs, DAE.LESS(tp), DAE.RCONST(upper), -1, NONE());
2114 then (str, u);
2115 end match;
2116
2117 196 algo := DAE.ALGORITHM_STMTS({DAE.STMT_ASSERT(cond, DAE.SCONST(msg), DAE.ASSERTIONLEVEL_ERROR, DAE.emptyElementSource)});
2118 98 assEq := BackendDAE.ALGORITHM(0, algo, DAE.emptyElementSource, DAE.EXPAND(), BackendDAE.EQ_ATTR_DEFAULT_UNKNOWN);
2119 end makeDomainAssert;
2120
2121 protected function makeInitialGuess
2122 input DAE.Type tp;
2123 input DAE.Exp iExp1;
2124 input DAE.Exp iExp2;
2125 output DAE.Exp oExp;
2126 protected
2127 DAE.Exp con, e;
2128 algorithm
2129 259 con := Expression.makePureBuiltinCall("initial", {}, tp);
2130 259 e := Expression.traverseExpBottomUp(iExp2, makeInitialGuess2, (iExp1, "pre", tp, true));
2131 259 oExp := Expression.traverseExpBottomUp(iExp2, makeInitialGuess2, (iExp1, "pre", tp, false));
2132 259 oExp := DAE.IFEXP(con, e, oExp);
2133 end makeInitialGuess;
2134
2135 protected function makeInitialGuess2
2136 input DAE.Exp iExp;
2137 input tuple<DAE.Exp, String, DAE.Type, Boolean> itpl;
2138 output DAE.Exp oExp;
2139 output tuple<DAE.Exp, String, DAE.Type, Boolean> otpl = itpl;
2140 algorithm
2141 oExp := match(iExp, itpl)
2142 local
2143 DAE.ComponentRef cr1,cr2;
2144 DAE.Type tp;
2145 String fun;
2146 DAE.Exp e;
2147
2148 case (DAE.CREF(componentRef=cr1), (DAE.CREF(componentRef=cr2), fun, tp, _))
2149 guard(ComponentReferenceBasics.crefEqual(cr1, cr2)) algorithm
2150 518 e := Expression.makePureBuiltinCall(fun, {iExp}, tp);
2151 then e;
2152
2153 case (_, (_, _, tp, true)) algorithm
2154 try
2155
3/4
✗ Branch 1 not taken.
✓ Branch 2 taken 132 times.
✓ Branch 3 taken 84 times.
✓ Branch 4 taken 48 times.
 132 SOME(e) := makeInitialGuess3(iExp, tp);
2156 else
2157 e := iExp;
2158 end try;
2159 then e;
2160
2161 else iExp;
2162 end match;
2163 end makeInitialGuess2;
2164
2165 protected function makeInitialGuess3
2166 input DAE.Exp iExp;
2167 input DAE.Type tp;
2168 output Option<DAE.Exp> oExp;
2169 algorithm
2170 oExp := match iExp
2171 local DAE.Exp e, con, o;
2172
2173 case DAE.CALL(path = Absyn.IDENT(name = "log"), expLst={e})
2174 algorithm
2175 12 con := DAE.RELATION(e, DAE.LESSEQ(tp), DAE.RCONST(0.0), -1, NONE());
2176 12 o := DAE.IFEXP(con, DAE.RCONST(-1/0.000000001), iExp);
2177 then SOME(o);
2178
2179 case DAE.CALL(path = Absyn.IDENT(name = "log10"), expLst={e})
2180 algorithm
2181 ✗ con := DAE.RELATION(e, DAE.LESSEQ(tp), DAE.RCONST(0.0), -1, NONE());
2182 ✗ o := DAE.IFEXP(con, DAE.RCONST(-1/0.000000001), iExp);
2183 then SOME(o);
2184
2185 case DAE.CALL(path = Absyn.IDENT(name = "sqrt"), expLst={e})
2186 algorithm
2187 ✗ con := DAE.RELATION(e, DAE.LESSEQ(tp), DAE.RCONST(0.0), -1, NONE());
2188 ✗ o := DAE.IFEXP(con, DAE.RCONST(0.0), iExp);
2189 then SOME(o);
2190
2191 case DAE.BINARY(exp2=e)
2192 algorithm
2193 36 con := DAE.RELATION(e, DAE.EQUAL(tp), DAE.RCONST(0.0), -1, NONE());
2194 36 o := DAE.IFEXP(con, DAE.RCONST(1.0), iExp);
2195 then SOME(o);
2196
2197 else NONE();
2198
2199 end match;
2200
2201 end makeInitialGuess3;
2202
2203 protected function helpInvCos
2204 input DAE.Exp acosy;
2205 input DAE.Exp x;
2206 input DAE.Type tp;
2207 input Boolean neg;
2208 output DAE.Exp k;
2209 algorithm
2210
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 33 times.
 66 k := if neg then
2211 Expression.expAdd(x,acosy)
2212 else
2213 Expression.expSub(x,acosy);
2214 66 k := Expression.makeDiv(k, Expression.expMul(DAE.RCONST(2.0), DAE.PI));
2215 66 k := Expression.makePureBuiltinCall("$_round",{k},tp);
2216
2217 end helpInvCos;
2218
2219 protected function helpInvSin
2220 input DAE.Exp asiny;
2221 input DAE.Exp x;
2222 input DAE.Type tp;
2223 input Boolean neg;
2224 output DAE.Exp k;
2225 algorithm
2226
2/2
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 27 times.
 54 k := if neg then
2227 Expression.expAdd(x,asiny)
2228 else
2229 Expression.expSub(x,asiny);
2230 54 k := Expression.makeDiv(k, Expression.expMul(DAE.RCONST(2.0), DAE.PI));
2231
2/2
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 27 times.
 54 if neg then
2232 27 k := Expression.expSub(k, DAE.RCONST(0.5));
2233 end if;
2234 54 k := Expression.makePureBuiltinCall("$_round",{k},tp);
2235 end helpInvSin;
2236
2237 protected function helpInvCos2
2238 input DAE.Exp k;
2239 input DAE.Exp acosy;
2240 input DAE.Type tp;
2241 input Boolean neg;
2242 output DAE.Exp x;
2243 algorithm
2244
2245
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 33 times.
 66 x := if neg then Expression.negate(acosy) else acosy;
2246 66 x := Expression.expAdd(x, Expression.expMul(k, Expression.expMul(DAE.RCONST(2.0), DAE.PI)));
2247
2248 end helpInvCos2;
2249
2250 protected function helpInvSin2
2251 input DAE.Exp k;
2252 input DAE.Exp asiny;
2253 input DAE.Type tp;
2254 input Boolean neg;
2255 output DAE.Exp x;
2256 protected
2257 DAE.Exp e;
2258 algorithm
2259
2260
2/2
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 27 times.
 54 x := if neg then Expression.negate(asiny) else asiny;
2261 54 e := Expression.expMul(k, Expression.expMul(DAE.RCONST(2.0), DAE.PI));
2262
2/2
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 27 times.
 54 e := if neg then Expression.expAdd(e, DAE.PI) else e;
2263 54 x := Expression.expAdd(x, e);
2264
2265 end helpInvSin2;
2266
2267 protected function helpInvCos3
2268 input DAE.Exp x1;
2269 input DAE.Exp x2;
2270 input DAE.Exp x;
2271 input DAE.Type tp;
2272 output DAE.Exp y;
2273 protected
2274 DAE.Exp diffx1 = absDiff(x1,x,tp);
2275 DAE.Exp diffx2 = absDiff(x2,x,tp);
2276 DAE.Exp con = DAE.RELATION(diffx1, DAE.LESS(tp), diffx2, -1, NONE());
2277 algorithm
2278 116 con := Expression.makeNoEvent(con);
2279 116 y := DAE.IFEXP(con, x1, x2);
2280 end helpInvCos3;
2281
2282 protected function absDiff
2283 input DAE.Exp x;
2284 input DAE.Exp y;
2285 input DAE.Type tp;
2286 output DAE.Exp z;
2287 algorithm
2288 232 z := Expression.expSub(x,y);
2289 232 z := Expression.makePureBuiltinCall("abs",{z},tp);
2290 end absDiff;
2291
2292 annotation(__OpenModelica_Interface="backend");
2293 end ExpressionSolve;
2294