Automatic state selection improved

- Models with lambda variables allowed; lambda variables are ignored when selecting non-states
- _x and _derx vectors of StateSelection printed, when an error occurs.

Author: MartinOtter

src/symbolic/DAEquations/BasicStructuralTransform.jl

@@ -994,7 +994,7 @@ function analyzeStructurally(equations, params, unknowns_indices, deriv, unknown
     for i in 1:length(IG)-length(Gsolvable)
         push!(Gsolvable, [])  # This is conservative since highest derivative is not marked as solvable.
     end
-    
+
     if automaticStateSelection 
 		startValues = []
 		fixedFlags = []
@@ -1004,7 +1004,7 @@ function analyzeStructurally(equations, params, unknowns_indices, deriv, unknown
 		end
 		@show startValues
 		@show fixedFlags
-	
+	    println(... 2)
         components = Array{Array{Int64,1},1}(BLT(IG, assign))
         vNames = makeList(unknownsNames, 1:length(assign), Avar) # ::Vector{String}
         eqGraph = StateSelection.getSortedEquationGraph(IG, Gsolvable, components, assign, Avar, Bequ, vNames, withStabilization=false)
@@ -1016,17 +1016,17 @@ function analyzeStructurally(equations, params, unknowns_indices, deriv, unknown
         newAssign = copy(assign)
         # newG = [newG[i] for i in 1:length(newG) if Bequ[i] == 0]
     end
-  
+
     vActive = fill(true, length(Avar))
     if ! automaticStateSelection
         vActive[statesIndices] .= false
     else
-		stateIndices = eqGraph.Vx
+		stateIndices = eqGraph.Vx[ findall(x->x>0, eqGraph.Vx) ]   # Only use Vx indices that are > 0 (= no lambda variables)
 		@show names Avar stateIndices
 		@show Gsolvable
 		for i in 1:length(stateIndices)
 			j = stateIndices[i]
-			if j in Avar # selected state is a derivative
+			if j in Avar  # selected state is a derivative, but not a lambda variable and not an algebraic variable
 				# search for corresponding variable
 				alias = 0
 				for ie in 1:length(G)
@@ -1049,7 +1049,7 @@ function analyzeStructurally(equations, params, unknowns_indices, deriv, unknown
         realStates = [GetField(This(), Symbol(name)) for name in stateNames]
         setRealStates(realStates)
     end
-    
+
     if log
         println("\nNot active:")
         for i in 1:length(vActive)

src/symbolic/StateSelection.jl

@@ -474,9 +474,12 @@ Return the sorted equation graph with selection of states and dummy states.
 function getSortedEquationGraph(G, Gsolvable, BLT, assign, A, B, VNames; withStabilization::Bool=true)
     eqInit = StateSelection.SortedEquationGraph(G, BLT, assign, A, B, VNames)
     eq = getSortedEquationGraph!(eqInit, Gsolvable)
+    #printSortedEquationGraph(eq; equations=false)
+    #println("\n")
 
     if !withStabilization
         # Check whether lambda and/or mue variables are present
+#=
         lambdaVariables = String[]
         for i in eachindex(eq.Vx)
             vx = eq.Vx[i]
@@ -489,25 +492,41 @@ function getSortedEquationGraph(G, Gsolvable, BLT, assign, A, B, VNames; withSta
 
         if length(lambdaVariables) > 0 || eq.nmue > 0
             # Code generator not prepared to handle DAE stabilization -> Trigger error
+            printSortedEquationGraph(eq; equations=false)
+
             constraintVariables = String[]             
             for i in eq.ider1n1
                 vc = eq.Vx[ eq.VxRev[ i ] ]
                 push!(constraintVariables, eq.VNames[ vc ])
             end
 
             if length(lambdaVariables) > 0
-                error("... Automatic state selection is not possible, because the code generator does not\n",
+                error("\n... Automatic state selection is not possible, because the code generator does not\n",
                       "    yet support DAE stabilization of higher index systems.\n",
                       "    Number of lambda variables: ", length(lambdaVariables), ", number of mue variables: ", eq.nmue, "\n",
-                      "    The following (lambda) variables are the reason: ", lambdaVariables, "\n",
-                      "    The following (state) variables are the reason: ", constraintVariables, ".")
+                      "    The following (lambda) variables are the reason: ", lambdaVariables, ".")
             elseif eq.nmue > 0
-                error("... Automatic state selection is not possible, because the code generator does not\n",
+                error("\n... Automatic state selection is not possible, because the code generator does not\n",
                       "    yet support the generation of stabilizing equations.\n",
                       "    Number of lambda variables: ", length(lambdaVariables), ", number of mue variables: ", eq.nmue, "\n",
                       "    The following (state) variables are the reason: ", constraintVariables, ".")      
             end
         end
+=#
+        if eq.nmue > 0
+            # Code generator not prepared to handle DAE stabilization -> Trigger error
+            printSortedEquationGraph(eq; equations=false)
+
+            constraintVariables = String[]             
+            for i in eq.ider1n1
+                vc = eq.Vx[ eq.VxRev[ i ] ]
+                push!(constraintVariables, eq.VNames[ vc ])
+            end
+            error("\n... Automatic state selection is not possible, because the code generator does not\n",
+                  "    yet support the generation of stabilizing equations.\n",
+                  "    Number of mue variables: ", eq.nmue, "\n",
+                  "    The following (state) variables are the reason: ", constraintVariables, ".")      
+        end
 
         # Determine differentiated variables in the original equations that are no states (dummy states)
     end
@@ -700,11 +719,12 @@ end
 
 
 """
-    printSortedEquationGraph(eq)
+    printSortedEquationGraph(eq; equations=true)
 
-Print information about the sorted equation graph
+Print information about the sorted equation graph.
+If `equations = false`, the equations are not printed.
 """
-function printSortedEquationGraph(eq::SortedEquationGraph)  
+function printSortedEquationGraph(eq::SortedEquationGraph; equations::Bool=true)  
    # Determine original number of equations (that are not differentiated)
     neqOrig = 0
     i = 0
@@ -786,6 +806,10 @@ function printSortedEquationGraph(eq::SortedEquationGraph)
         print("\n")
     end
 
+    if !equations
+        return
+    end
+
     # Print sorted equations
     println("\n  Sorted equations (length(_r) = ", length(eq.Er), ", nc = ", eq.nc, "):")