Skip to content
This repository was archived by the owner on May 18, 2019. It is now read-only.

Commit 275ccba

Browse files
ptaeuberOpenModelica-Hudson
ptaeuber
authored and
OpenModelica-Hudson
committed
Introduce new global homotopy approach
The new approach finds the smallest homotopy iteration loop in the system and solves it with an adaptive lambda step size. Activate with: --initOptModules+=generateHomotopyComponents (in combination with: --homotopyApproach=global) see ticket:2266
1 parent 62dc75d commit 275ccba

File tree

17 files changed

+1217
-557
lines changed

17 files changed

+1217
-557
lines changed

‎Compiler/BackEnd/BackendDAEOptimize.mo‎

Lines changed: 293 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5816,7 +5816,7 @@ algorithm
58165816
end for;
58175817
end inlineHomotopy;
58185818

5819-
public function inlineHomotopy2
5819+
protected function inlineHomotopy2
58205820
input BackendDAE.Equation inEq;
58215821
input Boolean inFoundHomotopy;
58225822
output BackendDAE.Equation outEq;
@@ -5855,5 +5855,297 @@ algorithm
58555855
end match;
58565856
end replaceHomotopyWithLambdaExpression;
58575857

5858+
// =============================================================================
5859+
// section for initOptModule >>generateHomotopyComponents<<
5860+
//
5861+
// =============================================================================
5862+
5863+
public function generateHomotopyComponents " finds the smallest homotopy loop and creates a component out of it
5864+
author: ptaeuber 2017"
5865+
input BackendDAE.BackendDAE inDAE;
5866+
output BackendDAE.BackendDAE outDAE = inDAE;
5867+
protected
5868+
BackendDAE.StrongComponents comps;
5869+
BackendDAE.EqSystems newEqSystems={};
5870+
array<Integer> ass1, ass2;
5871+
algorithm
5872+
for syst in outDAE.eqs loop
5873+
BackendDAE.MATCHING(ass1=ass1, ass2=ass2, comps=comps) := syst.matching;
5874+
(comps, syst) := traverseStrongComponentsForHomotopyLoop(comps, syst);
5875+
syst.matching := BackendDAE.MATCHING(ass1=ass1, ass2=ass2, comps=comps);
5876+
newEqSystems := syst::newEqSystems;
5877+
end for;
5878+
outDAE.eqs := listReverse(newEqSystems);
5879+
end generateHomotopyComponents;
5880+
5881+
protected function traverseStrongComponentsForHomotopyLoop " traverses all the strong components and finds the smallest homotopy loop
5882+
author: ptaeuber 2017"
5883+
input output BackendDAE.StrongComponents comps;
5884+
input output BackendDAE.EqSystem system;
5885+
protected
5886+
Integer nComps, compIndex=0, homotopyLoopBeginning=0, homotopyLoopEnd=0;
5887+
BackendDAE.StrongComponents preHomotopyComponents, homotopyComponents, postHomotopyComponents;
5888+
BackendDAE.StrongComponent homotopyComponent;
5889+
BackendDAE.Var lambda;
5890+
Integer lambdaIdx;
5891+
algorithm
5892+
nComps := listLength(comps);
5893+
for comp in comps loop
5894+
compIndex := compIndex + 1;
5895+
_ := match(comp)
5896+
local
5897+
Integer eqnIndex, varIndex;
5898+
list<Integer> eqnIndexes, varIndexes, resEqnIndexes, tVarIndexes, innerEqnIndexes;
5899+
list<list<Integer>> innerVarIndexesLst;
5900+
BackendDAE.InnerEquations innerEquations;
5901+
BackendDAE.Equation eqn;
5902+
list<BackendDAE.Equation> eqnLst;
5903+
Boolean hasHomotopy;
5904+
5905+
case(BackendDAE.SINGLEEQUATION(eqn=eqnIndex, var=varIndex))
5906+
equation
5907+
eqn = BackendEquation.get(system.orderedEqs, eqnIndex);
5908+
(_, hasHomotopy) = BackendEquation.traverseExpsOfEquation(eqn, BackendDAEUtil.containsHomotopyCall, false);
5909+
if hasHomotopy then
5910+
homotopyLoopEnd = compIndex;
5911+
if (homotopyLoopBeginning == 0) then
5912+
homotopyLoopBeginning = compIndex;
5913+
end if;
5914+
end if;
5915+
then();
5916+
5917+
case(BackendDAE.EQUATIONSYSTEM(eqns=eqnIndexes, vars=varIndexes))
5918+
equation
5919+
if (homotopyLoopBeginning == 0) then
5920+
eqnLst = BackendEquation.getList(eqnIndexes, system.orderedEqs);
5921+
(_, hasHomotopy) = BackendEquation.traverseExpsOfEquationList(eqnLst, BackendDAEUtil.containsHomotopyCall, false);
5922+
5923+
if hasHomotopy then
5924+
homotopyLoopBeginning = compIndex;
5925+
homotopyLoopEnd = compIndex;
5926+
end if;
5927+
else
5928+
homotopyLoopEnd = compIndex;
5929+
end if;
5930+
then();
5931+
5932+
case(BackendDAE.SINGLEARRAY(eqn=eqnIndex, vars=varIndexes))
5933+
equation
5934+
eqn = BackendEquation.get(system.orderedEqs, eqnIndex);
5935+
(_, hasHomotopy) = BackendEquation.traverseExpsOfEquation(eqn, BackendDAEUtil.containsHomotopyCall, false);
5936+
5937+
if hasHomotopy then
5938+
homotopyLoopEnd = compIndex;
5939+
if (homotopyLoopBeginning == 0) then
5940+
homotopyLoopBeginning = compIndex;
5941+
end if;
5942+
end if;
5943+
then();
5944+
5945+
case(BackendDAE.SINGLEALGORITHM(eqn=eqnIndex, vars=varIndexes))
5946+
equation
5947+
eqn = BackendEquation.get(system.orderedEqs, eqnIndex);
5948+
(_, hasHomotopy) = BackendEquation.traverseExpsOfEquation(eqn, BackendDAEUtil.containsHomotopyCall, false);
5949+
5950+
if hasHomotopy then
5951+
homotopyLoopEnd = compIndex;
5952+
if (homotopyLoopBeginning == 0) then
5953+
homotopyLoopBeginning = compIndex;
5954+
end if;
5955+
end if;
5956+
then();
5957+
5958+
case(BackendDAE.SINGLECOMPLEXEQUATION(eqn=eqnIndex, vars=varIndexes))
5959+
equation
5960+
eqn = BackendEquation.get(system.orderedEqs, eqnIndex);
5961+
(_, hasHomotopy) = BackendEquation.traverseExpsOfEquation(eqn, BackendDAEUtil.containsHomotopyCall, false);
5962+
5963+
if hasHomotopy then
5964+
homotopyLoopEnd = compIndex;
5965+
if (homotopyLoopBeginning == 0) then
5966+
homotopyLoopBeginning = compIndex;
5967+
end if;
5968+
end if;
5969+
then();
5970+
5971+
case(BackendDAE.SINGLEWHENEQUATION(eqn=eqnIndex, vars=varIndexes))
5972+
equation
5973+
eqn = BackendEquation.get(system.orderedEqs, eqnIndex);
5974+
(_, hasHomotopy) = BackendEquation.traverseExpsOfEquation(eqn, BackendDAEUtil.containsHomotopyCall, false);
5975+
5976+
if hasHomotopy then
5977+
homotopyLoopEnd = compIndex;
5978+
if (homotopyLoopBeginning == 0) then
5979+
homotopyLoopBeginning = compIndex;
5980+
end if;
5981+
end if;
5982+
then();
5983+
5984+
case(BackendDAE.SINGLEIFEQUATION(eqn=eqnIndex, vars=varIndexes))
5985+
equation
5986+
eqn = BackendEquation.get(system.orderedEqs, eqnIndex);
5987+
(_, hasHomotopy) = BackendEquation.traverseExpsOfEquation(eqn, BackendDAEUtil.containsHomotopyCall, false);
5988+
5989+
if hasHomotopy then
5990+
homotopyLoopEnd = compIndex;
5991+
if (homotopyLoopBeginning == 0) then
5992+
homotopyLoopBeginning = compIndex;
5993+
end if;
5994+
end if;
5995+
then();
5996+
5997+
case(BackendDAE.TORNSYSTEM(strictTearingSet=BackendDAE.TEARINGSET(residualequations=resEqnIndexes, tearingvars=tVarIndexes, innerEquations=innerEquations)))
5998+
equation
5999+
if (homotopyLoopBeginning == 0) then
6000+
eqnLst = BackendEquation.getList(resEqnIndexes, system.orderedEqs);
6001+
(_, hasHomotopy) = BackendEquation.traverseExpsOfEquationList(eqnLst, BackendDAEUtil.containsHomotopyCall, false);
6002+
6003+
if not hasHomotopy then
6004+
(innerEqnIndexes, innerVarIndexesLst,_) = List.map_3(innerEquations, BackendDAEUtil.getEqnAndVarsFromInnerEquation);
6005+
eqnLst = BackendEquation.getList(innerEqnIndexes, system.orderedEqs);
6006+
(_, hasHomotopy) = BackendEquation.traverseExpsOfEquationList(eqnLst, BackendDAEUtil.containsHomotopyCall, false);
6007+
end if;
6008+
6009+
if hasHomotopy then
6010+
homotopyLoopBeginning = compIndex;
6011+
homotopyLoopEnd = compIndex;
6012+
end if;
6013+
else
6014+
homotopyLoopEnd = compIndex;
6015+
end if;
6016+
then();
6017+
end match;
6018+
end for;
6019+
6020+
if homotopyLoopBeginning > 0 then
6021+
// Add homotopy lambda to system
6022+
lambda := BackendDAE.VAR(ComponentReference.makeCrefIdent(BackendDAE.homotopyLambda, DAE.T_REAL_DEFAULT, {}), BackendDAE.VARIABLE(), DAE.BIDIR(), DAE.NON_PARALLEL(), DAE.T_REAL_DEFAULT, NONE(), NONE(), {}, DAE.emptyElementSource, NONE(), NONE(), DAE.BCONST(false), NONE(), DAE.NON_CONNECTOR(), DAE.NOT_INNER_OUTER(), true);
6023+
system.orderedVars := BackendVariable.addVar(lambda, system.orderedVars);
6024+
lambdaIdx := BackendVariable.varsSize(system.orderedVars);
6025+
6026+
(preHomotopyComponents, homotopyComponents, postHomotopyComponents) := getHomotopyComponents(List.intRange(nComps), comps, homotopyLoopBeginning, homotopyLoopEnd);
6027+
6028+
homotopyComponent := createOneHomotopyComponent(homotopyComponents, system, lambdaIdx);
6029+
6030+
comps := homotopyComponent::postHomotopyComponents;
6031+
comps := listAppend(preHomotopyComponents, comps);
6032+
end if;
6033+
end traverseStrongComponentsForHomotopyLoop;
6034+
6035+
6036+
protected function getHomotopyComponents " divides the components into pre-homotopy, homotopy and post-homotopy parts
6037+
author: ptaeuber 2017"
6038+
input list<Integer> componentIndexes;
6039+
input BackendDAE.StrongComponents components;
6040+
input Integer homotopyLoopBeginning;
6041+
input Integer homotopyLoopEnd;
6042+
input output BackendDAE.StrongComponents outPreHomotopyComponents = {};
6043+
input output BackendDAE.StrongComponents outHomotopyComponents = {};
6044+
input output BackendDAE.StrongComponents outPostHomotopyComponents = {};
6045+
algorithm
6046+
(outPreHomotopyComponents, outHomotopyComponents, outPostHomotopyComponents) := match(componentIndexes, components)
6047+
local
6048+
Integer i;
6049+
list<Integer> indexes;
6050+
BackendDAE.StrongComponent comp;
6051+
BackendDAE.StrongComponents comps;
6052+
6053+
case(i::{}, comp::{})
6054+
equation
6055+
if intLt(i, homotopyLoopBeginning) then
6056+
outPreHomotopyComponents = comp::outPreHomotopyComponents;
6057+
elseif intGt(i, homotopyLoopEnd) then
6058+
outPostHomotopyComponents = comp::outPostHomotopyComponents;
6059+
else
6060+
outHomotopyComponents = comp::outHomotopyComponents;
6061+
end if;
6062+
then (listReverse(outPreHomotopyComponents), listReverse(outHomotopyComponents), listReverse(outPostHomotopyComponents));
6063+
6064+
case(i::indexes, comp::comps)
6065+
equation
6066+
if intLt(i, homotopyLoopBeginning) then
6067+
outPreHomotopyComponents = comp::outPreHomotopyComponents;
6068+
elseif intGt(i, homotopyLoopEnd) then
6069+
outPostHomotopyComponents = comp::outPostHomotopyComponents;
6070+
else
6071+
outHomotopyComponents = comp::outHomotopyComponents;
6072+
end if;
6073+
then getHomotopyComponents(indexes, comps, homotopyLoopBeginning, homotopyLoopEnd, outPreHomotopyComponents, outHomotopyComponents, outPostHomotopyComponents);
6074+
end match;
6075+
end getHomotopyComponents;
6076+
6077+
protected function createOneHomotopyComponent " creates one BackendDAE.TORNSYSTEM out of all homotopy components
6078+
author: ptaeuber"
6079+
input BackendDAE.StrongComponents homotopyComponents;
6080+
input BackendDAE.EqSystem inSystem;
6081+
input Integer lambdaIdx;
6082+
output BackendDAE.StrongComponent outHomotopyComponent;
6083+
protected
6084+
BackendDAE.InnerEquations newInnerEquations = {};
6085+
list<Integer> newResEquations = {};
6086+
list<Integer> newIterationVars = {};
6087+
Boolean isMixed = false;
6088+
algorithm
6089+
for comp in homotopyComponents loop
6090+
(newInnerEquations, newResEquations, newIterationVars) := match(comp)
6091+
local
6092+
Integer eqnIndex, varIndex;
6093+
list<Integer> eqnIndexes, varIndexes, resEqnIndexes, tVarIndexes;
6094+
BackendDAE.InnerEquations innerEquations;
6095+
BackendDAE.InnerEquation newInnerEquation;
6096+
Boolean mixedSystem;
6097+
6098+
case(BackendDAE.SINGLEEQUATION(eqn=eqnIndex, var=varIndex))
6099+
equation
6100+
newInnerEquation = BackendDAE.INNEREQUATION(eqn=eqnIndex, vars={varIndex});
6101+
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
6102+
6103+
case(BackendDAE.EQUATIONSYSTEM(eqns=eqnIndexes, vars=varIndexes, mixedSystem=mixedSystem))
6104+
equation
6105+
if mixedSystem then
6106+
isMixed = true;
6107+
end if;
6108+
then (newInnerEquations, listAppend(newResEquations, eqnIndexes), listAppend(newIterationVars, varIndexes));
6109+
6110+
case(BackendDAE.SINGLEARRAY(eqn=eqnIndex, vars=varIndexes))
6111+
equation
6112+
newInnerEquation = BackendDAE.INNEREQUATION(eqn=eqnIndex, vars=varIndexes);
6113+
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
6114+
6115+
case(BackendDAE.SINGLEALGORITHM(eqn=eqnIndex, vars=varIndexes))
6116+
equation
6117+
newInnerEquation = BackendDAE.INNEREQUATION(eqn=eqnIndex, vars=varIndexes);
6118+
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
6119+
6120+
case(BackendDAE.SINGLECOMPLEXEQUATION(eqn=eqnIndex, vars=varIndexes))
6121+
equation
6122+
newInnerEquation = BackendDAE.INNEREQUATION(eqn=eqnIndex, vars=varIndexes);
6123+
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
6124+
6125+
case(BackendDAE.SINGLEWHENEQUATION(eqn=eqnIndex, vars=varIndexes))
6126+
equation
6127+
newInnerEquation = BackendDAE.INNEREQUATION(eqn=eqnIndex, vars=varIndexes);
6128+
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
6129+
6130+
case(BackendDAE.SINGLEIFEQUATION(eqn=eqnIndex, vars=varIndexes))
6131+
equation
6132+
newInnerEquation = BackendDAE.INNEREQUATION(eqn=eqnIndex, vars=varIndexes);
6133+
then (newInnerEquation::newInnerEquations, newResEquations, newIterationVars);
6134+
6135+
case(BackendDAE.TORNSYSTEM(strictTearingSet=BackendDAE.TEARINGSET(residualequations=resEqnIndexes, tearingvars=tVarIndexes, innerEquations=innerEquations), mixedSystem=mixedSystem))
6136+
algorithm
6137+
if mixedSystem then
6138+
isMixed := true;
6139+
end if;
6140+
for innerEquation in innerEquations loop
6141+
newInnerEquations := innerEquation::newInnerEquations;
6142+
end for;
6143+
then (newInnerEquations, listAppend(newResEquations, resEqnIndexes), listAppend(newIterationVars, tVarIndexes));
6144+
end match;
6145+
end for;
6146+
6147+
outHomotopyComponent := BackendDAE.TORNSYSTEM(BackendDAE.TEARINGSET(listAppend(newIterationVars,{lambdaIdx}), newResEquations, listReverse(newInnerEquations), BackendDAE.EMPTY_JACOBIAN()), NONE(), false, isMixed);
6148+
end createOneHomotopyComponent;
6149+
58586150
annotation(__OpenModelica_Interface="backend");
58596151
end BackendDAEOptimize;

0 commit comments

Comments
 (0)