Documentation of StateSelection.jl and associated tests improved + new test model for automatic state selection added

Author: MartinOtter

src/symbolic/StateSelection.jl

@@ -156,6 +156,9 @@ mutable struct SortedEquationGraph
     Gunknowns::Vector{Vector{Int}}  # Gunknowns[i] is the vector of variables of equation i that are treated as unknowns
     Gsolvable::Vector{Any}          # Gsolvable[i] contains the variables that can be
                                     # solved for equation i (used for tearing)
+                                    # Gsolvable is only needed for the non-differentiated equations 
+                                    # (tearing variables of differentiated equations are automatically decuded from 
+                                    # the non-differentiated equations)
     eConstraintsVec::Vector{Vector{Vector{Int}}}   # ec = eConstraintsVec[i]: constraint sets of BLT block i; 
                                                    #      ec[1]  : lowest derivative equations
                                                    #      ec[end]: highest derivative equations (= BLT block i)
@@ -264,10 +267,10 @@ mutable struct SortedEquationGraph
             end
         end
 
-        if lowerDerivativeEquationsInBLT
-            println("\n... Warning from SortedEquationGraph(..) (in StateSelection.jl):",
-                "\n... BLT blocks with lower derivative equations have been ignored.\n")  
-        end 
+        #if lowerDerivativeEquationsInBLT
+        #    println("\n... Warning from SortedEquationGraph(..) (in StateSelection.jl):",
+        #        "\n... BLT blocks with lower derivative equations have been ignored.\n")  
+        #end 
 
       # Tearing information
         td = TraverseDAG(Gunknowns, length(VNames))
@@ -445,7 +448,8 @@ Return the sorted equation graph with selection of states and dummy states.
 
 # Input arguments
 - `Gorigin`: Gorigin[i] is the vector of variable indices of equation i.
-- `Gsolvable`: Gsolvable[i] is a subset of Gorigin[i] and contains the unknowns that can be explicitly solved for (used for tearing)
+- `Gsolvable`: Gsolvable[i] is a subset of Gorigin[i] and contains the unknowns that can be explicitly solved for (used for tearing).
+   Gsolvable is only needed for the non-differentiated equations (all differentiated equations j can have Gsolvable[j]=[])
 - `BLT`: BLT[i] is the vector of equations belonging to BLT-block i
 - `assign`: ei = assign[vi] is equation ei assigned to variable vi
 - `A`: A-Vector of Pantelides: `A[i] = if der(v[i]) == v[k] then k else 0` where `v[i]` is variable `i`. 
@@ -520,7 +524,8 @@ Constructs all information for the sorted equation graph.
 Input arguments:
 
 - `eqInit`: Initialized data structure instantiated with SortedEquationGraph()
-- `Gsolvable`: Subset of eqInit.Gunknowns that contains the unknowns that can be solved for (used for tearing)
+- `Gsolvable`: Subset of eqInit.Gunknowns that contains the unknowns that can be solved for (used for tearing).
+   Gsolvable is only needed for the non-differentiated equations (all differentiated equations j can have Gsolvable[j]=[])
 
 Output arguments in sortedEquationGraph:
 

test/models/TestAutomaticStateSelection.jl

@@ -83,8 +83,6 @@ result = simulate(ParallelCapacitors2, 1; logTranslation=true, logSimulation=log
 plot(result, ("C1.v", "C2.v"), figure=4)
 
 
-# Simulates with tearing = false and u2.state=false
-# Translation error if tearing = true and u2.state=false
 @model ParallelCapacitors2b begin
     C1 = 1e-3
     C2 = 2e-3
@@ -105,4 +103,26 @@ result = simulate(ParallelCapacitors2b , 1.0; logTranslation=true, logSimulation
 plot(result, ("u1", "u2", "i1", "v1"), figure=6)
 
 
+@model TwoInertiasConnectedViaIdealGearWithPositionConstraints begin
+    J1 = 2.0
+    J2 = 3.0
+    ratio = 4.0
+    phi1 = Float(size=(), start=1.0)
+    phi2 = Float(size=(), start=1.0)
+    w1 = Float(size=(), start=1.0)
+    w2 = Float(size=(), start=1.0)
+    tau = Float(size=())
+    t   = Float(size=())
+    @equations begin
+        der(t) = 1.0
+        w1 = der(phi1)
+        w2 = der(phi2)
+        J1*der(w1) = sin(t) - ratio*tau
+        J2*der(w2) = tau
+        phi1 = ratio*phi2
+    end
+end 
+result = simulate(TwoInertiasConnectedViaIdealGearWithPositionConstraints, 3.0; logTranslation=true, logSimulation=logSimulation, tearing=true, removeSingularities=true, automaticStateSelection=true)
+plot(result, [("phi1", "phi2"), ("w1", "w2")], figure=7)
+
 end

test/symbolic/TestStateSelectionAlgorithm.jl

@@ -148,7 +148,9 @@ vNames:
 
         # First version where all unknowns can be solved for
         println("\n\n... Test two coupled inertias (all unknowns can be solved for)")
-        Gsolvable = Gorigin
+        # Gsolvable = Gorigin
+        Gsolvable = StateSelection.newRaggedIntMatrix(length(Gorigin))
+        Gsolvable[1:6] = Gorigin[1:6]
         eq = getSortedEquationGraph(Gorigin, Gsolvable, BLT, assign, Avar, Bequ, VNames)
         if printSortedEquations
            printSortedEquationGraph(eq)
@@ -157,7 +159,7 @@ vNames:
 
 
         println("\n\n... Test two coupled inertias (only a subset of unknowns can be solved for)")
-        Gsolvable = Any[Any[3,7],Any[1],Any[6,9],Any[4],Any[2,5],Any[1,4],Any[7,9],[8,11],[10,12],[12,11]]
+        Gsolvable = Any[Any[3,7],Any[1],Any[6,9],Any[4],Any[2,5],Any[1,4],Any[],[],[],[]]
         eq = getSortedEquationGraph(Gorigin, Gsolvable, BLT, assign, Avar, Bequ, VNames)
         if printSortedEquations
            printSortedEquationGraph(eq)
@@ -214,7 +216,7 @@ vNames:
    9: (9)
 =#
 
-
+        println("\n\n... Test simple sliding mass model with Tearing")
         assign = [0,0,0,3,5,0,4,0,7,8,0,9]
         A = [6,8,10,0,0,7,0,9,0,0,12,0]
         B = [6,8,0,0,0,7,0,0,0]
@@ -235,6 +237,7 @@ vNames:
      [10,9],
      [12,11]]
 
+    #=
         Gsolvable = [[2],
              [3,8],
              [10,4,5],
@@ -244,6 +247,17 @@ vNames:
              [9],
              [10,9],
              [12]]
+    =#
+
+        Gsolvable = [[2],
+             [3,8],
+             [10,4,5],
+             [],
+             [5],
+             [],
+             [],
+             [],
+             [12]]
 
         eq = getSortedEquationGraph(G, Gsolvable, BLT, assign, A, B, VNames)
         if printSortedEquations
@@ -359,8 +373,9 @@ vNames:
 
 
         println("\n\n... Test Multi-Index DAE WITH tearing")
-        Gsolvable    = copy(G)
-        Gsolvable[8] = fill(0, 0)  # "0 = u8(t) + x8 - sin(x8)" cannot be solved for x8
+        # Gsolvable    = copy(G)
+        # Gsolvable[8] = fill(0, 0)  # "0 = u8(t) + x8 - sin(x8)" cannot be solved for x8
+        Gsolvable[1:7] = G[1:7]   
         eq = getSortedEquationGraph(G, Gsolvable, BLT, assign, A, B, VNames)
         if printSortedEquations
            printSortedEquationGraph(eq)