[BE] Add jacobian dependencies

  - Preparation for initial state selection, used for matching
  - The dependencies of a jacobian are all unknowns contained in it

Belonging to [master]:
  - OpenModelica/OpenModelica#136
  - #3048

Author: kabdelhak

Compiler/BackEnd/BackendDAE.mo

@@ -333,7 +333,7 @@ public constant EvaluationStages defaultEvalStages = EVALUATION_STAGES(false,fal
 
 public uniontype EquationAttributes
   record EQUATION_ATTRIBUTES
-    Boolean differentiated "true if the equation was differentiated, and should not differentiated again to avoid equal equations";
+    Boolean differentiated "true if the equation was differentiated, and should not be differentiated again to avoid equal equations";
     EquationKind kind;
     EvaluationStages evalStages;
   end EQUATION_ATTRIBUTES;
@@ -780,7 +780,8 @@ type SymbolicJacobian = tuple<BackendDAE,               // symbolic equation sys
                               String,                   // Matrix name
                               list<Var>,                // diff vars (independent vars)
                               list<Var>,                // diffed vars (residual vars)
-                              list<Var>                 // all diffed vars (residual vars + dependent vars)
+                              list<Var>,                // all diffed vars (residual vars + dependent vars)
+                              list< .DAE.ComponentRef>    // original dependent variables
                               >;
 
 public

Compiler/BackEnd/BackendDAEOptimize.mo

@@ -1463,7 +1463,7 @@ algorithm
         DAE.FunctionTree usedfuncs;
         BackendDAE.Shared shared;
 
-      case (SOME((bdae, _, _, _, _)), _, _)
+      case (SOME((bdae, _, _, _, _, _)), _, _)
         equation
           bdae = BackendDAEUtil.setFunctionTree(bdae, inFunctions);
           shared = bdae.shared;
@@ -1930,7 +1930,7 @@ algorithm
       /* TODO: implement/check for GENERIC_JACOBIAN */
       case BackendDAE.GENERIC_JACOBIAN(jacobian=NONE())
         then -1;
-      case BackendDAE.GENERIC_JACOBIAN(jacobian=SOME((_,_,vars1,vars2,_)))
+      case BackendDAE.GENERIC_JACOBIAN(jacobian=SOME((_,_,vars1,vars2,_,_)))
         guard
           listLength(vars1) == listLength(vars2)
         then listLength(vars1);

Compiler/BackEnd/BackendDAEUtil.mo

@@ -6267,9 +6267,9 @@ algorithm
         BackendDAE.BackendDAE bdae;
       case BackendDAE.EQUATIONSYSTEM(jac=BackendDAE.FULL_JACOBIAN(SOME(jac)))
         then traverseBackendDAEExpsJacobianEqn(jac, inFunc, arg);
-      case BackendDAE.EQUATIONSYSTEM(jac=BackendDAE.GENERIC_JACOBIAN(jacobian = SOME((bdae,_,_,_,_))))
+      case BackendDAE.EQUATIONSYSTEM(jac=BackendDAE.GENERIC_JACOBIAN(jacobian = SOME((bdae,_,_,_,_,_))))
         then traverseBackendDAEExps(bdae, inFunc, arg);
-      case BackendDAE.TORNSYSTEM(BackendDAE.TEARINGSET(jac=BackendDAE.GENERIC_JACOBIAN(jacobian = SOME((bdae,_,_,_,_)))))
+      case BackendDAE.TORNSYSTEM(BackendDAE.TEARINGSET(jac=BackendDAE.GENERIC_JACOBIAN(jacobian = SOME((bdae,_,_,_,_,_)))))
         then traverseBackendDAEExps(bdae, inFunc, arg);
       else arg;
     end match;
@@ -6325,7 +6325,7 @@ algorithm
       BackendDAE.BackendDAE bdae;
       Type_a arg;
     case ({}, _, _) then inTypeA;
-    case (BackendDAE.STATESET(jacobian = BackendDAE.GENERIC_JACOBIAN(jacobian = SOME((bdae,_,_,_,_))))::rest, _, _)
+    case (BackendDAE.STATESET(jacobian = BackendDAE.GENERIC_JACOBIAN(jacobian = SOME((bdae,_,_,_,_,_))))::rest, _, _)
       equation
         arg = traverseBackendDAEExps(bdae, inFunc, inTypeA);
       then
@@ -7052,6 +7052,9 @@ algorithm
     dae := BackendDAEOptimize.evaluateOutputsOnly(dae);
   end if;
 
+  //generate Jacobian for StateSets for initial state selection
+  dae := SymbolicJacobian.calculateStateSetsJacobians(dae);
+
   // generate system for initialization
   (outInitDAE, outInitDAE_lambda0_option, outRemovedInitialEquationLst, globalKnownVars) := Initialization.solveInitialSystem(dae);
   if Flags.isSet(Flags.WARN_NO_NOMINAL) then

Compiler/BackEnd/BackendDump.mo

@@ -2320,7 +2320,7 @@ algorithm
     then s;
   case(BackendDAE.GENERIC_JACOBIAN(jacobian=SOME(sJac),sparsePattern=sparsePattern))
     equation
-      ((dae,_,_,_,_)) = sJac;
+      ((dae,_,_,_,_, _)) = sJac;
       s = "GENERIC JACOBIAN:\n";
       dumpBackendDAE(dae,"Directional Derivatives System");
       dumpSparsityPattern(sparsePattern,"Sparse Pattern");
@@ -2351,7 +2351,7 @@ algorithm
       String s;
   case((SOME(sJac), sparsePattern, _))
     equation
-      ((dae,_,_,_,_)) = sJac;
+      ((dae,_,_,_,_, _)) = sJac;
       s = "GENERIC JACOBIAN:\n";
       dumpBackendDAE(dae,"Directional Derivatives System");
       dumpSparsityPattern(sparsePattern,"Sparse Pattern");

Compiler/BackEnd/BackendEquation.mo

@@ -609,6 +609,28 @@ algorithm
   (_, (_,cr_lst)) := traverseExpsOfEquation(e, Expression.traverseSubexpressionsHelper, (Expression.traversingComponentRefFinder, cr_lst));
 end traversingEquationCrefFinder;
 
+public function getCrefsFromEquations
+  input BackendDAE.EquationArray inEqns;
+  input BackendDAE.Variables inVars;
+  input BackendDAE.Variables inKnVars;
+  output list<DAE.ComponentRef> cr_lst;
+protected
+  HashTable.HashTable ht;
+  BackendDAE.Variables vars;
+  BackendDAE.Variables knownVars;
+algorithm
+  ht := HashTable.emptyHashTable();
+  (_, _, ht) := traverseEquationArray(inEqns, findUnknownCrefs, (inVars, inKnVars, ht));
+  cr_lst := BaseHashTable.hashTableKeyList(ht);
+end getCrefsFromEquations;
+
+protected function findUnknownCrefs
+  input output BackendDAE.Equation inEq;
+  input output tuple<BackendDAE.Variables, BackendDAE.Variables, HashTable.HashTable> extraArgs;
+algorithm
+  (_, ( _, extraArgs)) := traverseExpsOfEquation(inEq, Expression.traverseSubexpressionsHelper, (checkEquationsUnknownCrefsExp, extraArgs));
+end findUnknownCrefs;
+
 public function equationUnknownCrefs "author: Frenkel TUD 2012-05
   From the equation and a variable array return all
   variables in the equation an not in the variable array."
@@ -3182,6 +3204,5 @@ algorithm
   end match;
 end scalarComplexEquations;
 
-
 annotation(__OpenModelica_Interface="backend");
 end BackendEquation;

Compiler/BackEnd/FindZeroCrossings.mo

@@ -892,10 +892,11 @@ protected
   BackendDAE.BackendDAE jacBDAE;
   String name;
   list<BackendDAE.Var> seedVars, tmpVars, resultVars;
+  list<DAE.ComponentRef> depCrefs;
 algorithm
-  (jacBDAE, name, seedVars, tmpVars, resultVars) := symJac;
+  (jacBDAE, name, seedVars, tmpVars, resultVars, depCrefs) := symJac;
   jacBDAE := replaceZeroCrossingsJacBackend(jacBDAE, zeroCrossingLst, relationsLst, samplesLst, allVariables, globalKnownVars);
-  outSymJac := (jacBDAE, name, seedVars, tmpVars, resultVars);
+  outSymJac := (jacBDAE, name, seedVars, tmpVars, resultVars, depCrefs);
 end replaceZCExpinSymJacobian;
 
 protected function replaceZeroCrossingsJacBackend

Compiler/BackEnd/SymbolicJacobian.mo

@@ -1744,7 +1744,7 @@ try
     if Flags.isSet(Flags.JAC_DUMP2) then
       BackendDump.dumpSparsityPattern(sparsePattern, "FMI sparsity");
     end if;
-    outJacobianMatrixes := (SOME((emptyBDAE,"FMIDER",{},{},{})), sparsePattern, sparseColoring)::outJacobianMatrixes;
+    outJacobianMatrixes := (SOME((emptyBDAE,"FMIDER",{},{},{}, {})), sparsePattern, sparseColoring)::outJacobianMatrixes;
     outFunctionTree := inBackendDAE.shared.functionTree;
   else
     // prepare more needed variables
@@ -2001,7 +2001,7 @@ algorithm
     local
       BackendDAE.BackendDAE backendDAE, reducedDAE;
 
-      list<DAE.ComponentRef> comref_vars, comref_differentiatedVars;
+      list<DAE.ComponentRef> comref_vars, comref_differentiatedVars, dependencies;
 
       BackendDAE.Shared shared;
       BackendDAE.Variables  globalKnownVars, globalKnownVars1;
@@ -2043,9 +2043,10 @@ algorithm
         if Flags.isSet(Flags.JAC_DUMP2) then
           print("analytical Jacobians -> generated Jacobian DAE time: " + realString(clock()) + "\n");
         end if;
+        dependencies = calcJacobianDependencies((backendDAE, "", {}, {}, {}, {}));
 
      then
-        ((backendDAE, inName, inDiffVars, diffedVars, inVars), funcs);
+        ((backendDAE, inName, inDiffVars, diffedVars, inVars, dependencies), funcs);
     else
       equation
         Error.addInternalError("function createJacobian failed", sourceInfo());
@@ -2382,7 +2383,7 @@ algorithm
       BackendDAE.SymbolicJacobians rest;
       String name;
 
-    case (matrix as (SOME((_,name,_,_,_)), _, _))::_ guard
+    case (matrix as (SOME((_,name,_,_,_,_)), _, _))::_ guard
       stringEq(name, inJacobianName)
     then SOME(matrix);
 
@@ -2393,6 +2394,36 @@ algorithm
   end match;
 end getJacobianMatrixbyName;
 
+
+public function calcJacobianDependencies
+  input BackendDAE.SymbolicJacobian jacobian;
+  output list<DAE.ComponentRef> dependencies;
+protected
+  BackendDAE.EqSystem syst;
+  BackendDAE.Shared shared;
+algorithm
+  (BackendDAE.DAE({syst}, shared), _, _, _, _, _) := jacobian;
+  dependencies := BackendEquation.getCrefsFromEquations(syst.orderedEqs, syst.orderedVars, shared.globalKnownVars);
+end calcJacobianDependencies;
+
+public function getJacobianDependencies
+  input BackendDAE.Jacobian jacobian;
+  output list<DAE.ComponentRef> dependencies;
+algorithm
+  dependencies := match(jacobian)
+    case (BackendDAE.GENERIC_JACOBIAN(jacobian=SOME((_, _, _, _, _, dependencies))))
+    then dependencies;
+
+    case (BackendDAE.GENERIC_JACOBIAN(jacobian=NONE()))
+    then {};
+
+    else equation
+      Error.addInternalError("function getJacobianDependencies failed", sourceInfo());
+    then fail();
+
+  end match;
+end getJacobianDependencies;
+
 // =============================================================================
 // Module for to calculate strong component Jacobains
 //
@@ -2826,7 +2857,8 @@ algorithm
 
     // create known variables
     knvarLst1 := BackendEquation.equationsVars(eqns, globalKnownVars);
-    knvarLst2 := BackendEquation.equationsVars(eqns, inAllVars);
+    //knvarLst2 := BackendEquation.equationsVars(eqns, inAllVars);
+    knvarLst2 := {};
     // Create a list of known variables true *only* for this shared system
     globalKnownVars := BackendVariable.listVar2(knvarLst1,knvarLst2);
     // Remove inputs for the jacobian
@@ -3716,7 +3748,7 @@ protected
   BackendDAE.BackendDAE jacBDAE;
   String name;
 algorithm
-  (jacBDAE, name, _, _, _) := symbolicJacobian;
+  (jacBDAE, name, _, _, _, _) := symbolicJacobian;
   try
     _ := BackendDAEUtil.mapEqSystem(jacBDAE, checkForNonLinearStrongComponents_work);
     result := true;
@@ -3754,5 +3786,123 @@ algorithm
   end try;
 end checkForNonLinearStrongComponents_work;
 
+
+public function getFixedStatesForSelfdependentSets
+  " author: kabdelhak
+  Returns states to fix for initial problem in the case of selfdependent dynamic state sets"
+  input BackendDAE.StateSet stateSet;
+  input list<BackendDAE.Var> unfixedStates;
+  input Integer toFix;
+  output list<BackendDAE.Var> statesToFix;
+protected
+  list<tuple<Integer,BackendDAE.Var>> nonlinearCountLst = {};
+  Integer nonlinearCount;
+algorithm
+  _:= match(stateSet.jacobian)
+  local
+    BackendDAE.SymbolicJacobian sJac;
+    BackendDAE.BackendDAE dae;
+    list<BackendDAE.Var> diffVars;
+    String matrixName;
+  case (BackendDAE.GENERIC_JACOBIAN(jacobian=SOME(sJac))) algorithm
+    ((dae,matrixName,diffVars, _, _,_)) := sJac;
+    for var in unfixedStates loop
+      nonlinearCountLst := getNonlinearStateCount(var,unfixedStates,dae,matrixName)::nonlinearCountLst;
+    end for;
+  then 0;
+  end match;
+  statesToFix := fixedVarsFromNonlinearCount(nonlinearCountLst, toFix);
+end getFixedStatesForSelfdependentSets;
+
+protected function getNonlinearStateCount
+  input BackendDAE.Var state;
+  input list<BackendDAE.Var> diffVars;
+  input BackendDAE.BackendDAE dae;
+  input String matrixName;
+  output tuple<Integer,BackendDAE.Var> outTpl;
+protected
+algorithm
+  outTpl:=match(dae)
+  local
+    BackendDAE.EqSystems systs;
+    tuple<BackendDAE.Var,list<BackendDAE.Var>,Integer,String> tpl;
+    BackendDAE.Var outState;
+    Integer nonlinearCount = 0;
+  case BackendDAE.DAE(eqs=systs) algorithm
+    tpl := (state,diffVars,nonlinearCount,matrixName);
+    for syst in systs loop
+      _:= match(syst)
+      local
+        BackendDAE.EquationArray eqnarray;
+
+      case BackendDAE.EQSYSTEM(_,eqnarray,_,_,_,_,_,_) algorithm
+        tpl := BackendEquation.traverseEquationArray(eqnarray,getNonlinearStateCount0,tpl);
+      then 0;
+      end match;
+    end for;
+    (outState,_,nonlinearCount,_) := tpl;
+  then (nonlinearCount,outState);
+  end match;
+
+end getNonlinearStateCount;
+
+protected function getNonlinearStateCount0
+  input BackendDAE.Equation inEq;
+  input tuple<BackendDAE.Var,list<BackendDAE.Var>,Integer,String> inTpl;
+  output BackendDAE.Equation outEq;
+  output tuple<BackendDAE.Var,list<BackendDAE.Var>,Integer,String> outTpl;
+algorithm
+  outEq := inEq;
+  outTpl := match inEq
+  local
+    DAE.Exp exp, diffExp;
+    BackendDAE.Var state;
+    list<BackendDAE.Var> diffVars;
+    Integer nonlinearCount;
+    String matrixName;
+    DAE.ComponentRef seedVar;
+    list<DAE.Subscript> subs;
+  case BackendDAE.EQUATION(scalar=exp) algorithm
+    (state,diffVars,nonlinearCount,matrixName) := inTpl;
+    // Differentiate equation to look for nonlinear dependencies
+    seedVar := Differentiate.createSeedCrefName(BackendVariable.varCref(state),matrixName);
+    diffExp := Differentiate.differentiateExpSolve(exp,seedVar,NONE());
+    for var in diffVars loop
+      if not ComponentReference.crefEqual(var.varName, state.varName) and Expression.expContains(diffExp,Expression.crefExp(var.varName)) then
+        // Heuristic to punish fixed vars with a value of zero
+        if BackendVariable.varFixed(var) and Expression.isZero(BackendVariable.varStartValue(var)) then
+          nonlinearCount := nonlinearCount + 2;
+        else
+          nonlinearCount := nonlinearCount + 1;
+        end if;
+      end if;
+    end for;
+  then (state,diffVars,nonlinearCount,matrixName);
+  end match;
+end getNonlinearStateCount0;
+
+protected function fixedVarsFromNonlinearCount
+  input list<tuple<Integer,BackendDAE.Var>> tplLst;
+  input Integer toFix;
+  output list<BackendDAE.Var> fixedVars = {};
+protected
+  list<tuple<Integer,BackendDAE.Var>> sortedTplLst, strippedTplLst;
+  BackendDAE.Var fixVar;
+  Integer fixInt;
+algorithm
+  for tpl in tplLst loop
+     (fixInt,fixVar) := tpl;
+  end for;
+  // Sort by nonlinear count and take first N states to fix
+  sortedTplLst := List.sort(tplLst, Util.compareTupleIntGt);
+  strippedTplLst := List.firstN(sortedTplLst,toFix);
+  for tpl in strippedTplLst loop
+    (_,fixVar) := tpl;
+    fixVar.values := DAEUtil.setFixedAttr(fixVar.values,SOME(DAE.BCONST(true)));
+    fixedVars := fixVar::fixedVars;
+  end for;
+end fixedVarsFromNonlinearCount;
+
+
 annotation(__OpenModelica_Interface="backend");
 end SymbolicJacobian;