Handle automatic state selection when constrains on angles

Author: HildingElmqvist

src/Modia.jl

@@ -46,7 +46,7 @@ For more information, see (https://github.com/ModiaSim/Modia.jl/blob/master/READ
 module Modia
 
 const Version = "0.2.4-dev"
-const Date = "2019-02-13"
+const Date = "2019-02-24"
 
 
 #println(" \n\nWelcome to Modia - Dynamic MODeling and Simulation in julIA")

src/symbolic/DAEquations/BasicStructuralTransform.jl

@@ -996,9 +996,18 @@ function analyzeStructurally(equations, params, unknowns_indices, deriv, unknown
     end
 
     if automaticStateSelection 
+		startValues = []
+		fixedFlags = []
+		for (name, var) in unknowns
+			push!(startValues, var.start)
+			push!(fixedFlags, var.fixed)
+		end
+		@show startValues
+		@show fixedFlags
+	
         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)
+        eqGraph = StateSelection.getSortedEquationGraph(IG, Gsolvable, components, assign, Avar, Bequ, vNames, withStabilization=false)
     end
 
     if newStateSelection
@@ -1012,8 +1021,32 @@ function analyzeStructurally(equations, params, unknowns_indices, deriv, unknown
     if ! automaticStateSelection
         vActive[statesIndices] .= false
     else
-        vActive[eqGraph.Vx] .= false
-        realStates = [GetField(This(), name) for name in names[eqGraph.Vx]]
+		stateIndices = eqGraph.Vx
+		@show names Avar stateIndices
+		@show Gsolvable
+		for i in 1:length(stateIndices)
+			j = stateIndices[i]
+			if j in Avar # selected state is a derivative
+				# search for corresponding variable
+				alias = 0
+				for ie in 1:length(G)
+				    g = Gsolvable[ie]
+					if j in g && length(g) == 2 && sort(g) == sort(G[ie])
+						alias = if g[1] == j; g[2] else g[1] end
+						break;
+					end
+				end
+				@show alias i
+				stateIndices[i] = alias
+			end
+		end
+		
+		@show stateIndices
+        vActive[stateIndices] .= false
+        stateNames = [replace(replace(string(name), "(" => "_"), ")" => "") for name in names[stateIndices]]
+        @show stateNames
+		@show realStates
+        realStates = [GetField(This(), Symbol(name)) for name in stateNames]
         setRealStates(realStates)
     end
 

src/symbolic/DAEquations/SymbolicTransform.jl

@@ -589,22 +589,26 @@ function differentiate(e)
                 print("DIFFERENTIATE DERIVATIVE:    "); println(e)
         end
 
-        # @show e.base e.base.name
-        if true
+#        @show e e.base e.base.name realStates
+        state = GetField(This(), Symbol("der_" * string(e.base.name)))
+#        @show state        
+        i = findfirst(isequal(state), realStates)
+#        @show i
+        if i == notFound
             diff = GetField(This(), Symbol("der_der_" * string(e.base.name))) 
+#            @show diff
             if !noUnits 
                 diff = diff / Unitful.s
             end
             dummyDerivatives[Symbol("der_der_" * string(e.base.name))] = nothing
-        elseif e.base.name == Symbol("j.s")   # Hack!!!
-            diff = Der(GetField(This(), Symbol("j.v"))) # GetField(This(), Symbol("j.a")) #
         else
-            diff = Der(e)
+            diff = Der(state)
+            dummyDerivatives[Symbol("der_" * string(e.base.name))] = nothing
         end
 
         # @show e diff
         if logDifferentiateVariable
-                print("  Differentiated derivative: "); println(diff)
+            print("  Differentiated derivative: "); println(diff)
         end
 
     elseif e.head in [:(=), :(:=)] # Equation

test/models/TestAutomaticStateSelection.jl

@@ -16,8 +16,9 @@ else
     using Modia.Test
 end
 
-const logSimulation=false
+const logSimulation=true
 
+#=
 @model TwoConnectedInertias begin
     J1 = 2.0
     J2 = 3.0
@@ -101,17 +102,17 @@ plot(result, ("C1.v", "C2.v"), figure=4)
 end 
 result = simulate(ParallelCapacitors2b , 1.0; logTranslation=true, logSimulation=logSimulation, tearing=true, removeSingularities=true, automaticStateSelection=false)
 plot(result, ("u1", "u2", "i1", "v1"), figure=6)
+=#
 
-#=
 # Does not compile
 @model TwoInertiasConnectedViaIdealGearWithPositionConstraints begin
     J1 = 2.0
     J2 = 3.0
     ratio = 4.0
-    phi1 = Float(size=(), start=1.0)
+    phi1 = Float(size=(), start=1.0, fixed=true)
     phi2 = Float(size=(), start=1.0)
-    w1 = Float(size=(), start=1.0)
-    w2 = Float(size=(), start=1.0)
+    w1 = Float(size=(), start=1.0, fixed=true)
+    w2 = Float(size=(), start=1.0, fixed=true)
     tau = Float(size=())
     t   = Float(size=())
     @equations begin
@@ -123,9 +124,9 @@ plot(result, ("u1", "u2", "i1", "v1"), figure=6)
         phi1 = ratio*phi2
     end
 end 
-result = simulate(TwoInertiasConnectedViaIdealGearWithPositionConstraints, 3.0; logTranslation=true, logSimulation=logSimulation, tearing=true, removeSingularities=true, automaticStateSelection=true)
+result = simulate(TwoInertiasConnectedViaIdealGearWithPositionConstraints, 3.0; logTranslation=true, logSimulation=logSimulation, tearing=true, removeSingularities=false, automaticStateSelection=true)
 plot(result, [("phi1", "phi2"), ("w1", "w2")], figure=7)
-=#
+println("after plot")
 
 
 end