Skip to content

Commit

Permalink
Bug fixed, if algebraic equations appear differentiated + new optiona…
Browse files Browse the repository at this point in the history 
…l arguments var_is_state + logStates
  • Loading branch information
@MartinOtter
MartinOtter committed Feb 6, 2021
1 parent 7db826f commit 09a1c47
Showing 1 changed file with 63 additions and 20 deletions.
83 changes: 63 additions & 20 deletions src/StateSelection.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ const ERROR = 2
var_has_start = v_original -> false,
var_fixed = v_original -> nothing,
var_length = v_original -> 1,
var_is_state = v_original -> false,
equation = e -> ""
isSolvableEquation = (e_original,v_original) -> false,
isLinearEquation = (e_original,v_original) -> (false, false),
Expand Down Expand Up @@ -78,6 +79,9 @@ The functions need to have the following arguments:
- `var_length(v_original::Int)::Int`:\\
Return length of variable `v_original` from the original, undifferentiated equations.
- `var_is_state(v_original::Int)::Bool`:\\
Return true, if variable `v_original` is defined to be a state.
- `equation(e::Int)`:\\
Return equation as string or "", if equation is not provided.
Expand Down Expand Up @@ -113,6 +117,7 @@ struct StateSelectionFunctions
var_nominal::Function
var_unbounded::Function
var_length::Function
var_is_state::Function
equation::Function
isSolvableEquation::Function
isLinearEquation::Function
Expand All @@ -129,6 +134,7 @@ struct StateSelectionFunctions
var_nominal = v_original -> NaN,
var_unbounded = v_original -> false,
var_length = v_original -> 1,
var_is_state = v_original -> false,
equation = e -> "",
isSolvableEquation = (e_original,v_original) -> false,
isLinearEquation = (e_original,v_original) -> (false, false),
Expand All @@ -144,6 +150,7 @@ struct StateSelectionFunctions
var_nominal,
var_unbounded,
var_length,
var_is_state,
equation,
isSolvableEquation,
isLinearEquation,
Expand Down Expand Up @@ -207,7 +214,8 @@ showCodeWithoutComments(code) = println("code = ", replace(sprint(show,code), r"


"""
(eConstraints,vConstraints) = getConstraintSets(BLT_i, assignRev, Arev, Brev, showMessage)
(eConstraints,vConstraints) = getConstraintSets(BLT_i, assignRev, Arev, Brev,
var_is_state, showMessage)
Determines the set of equations/constraints and their unknowns that are associated with
the highest derivative level BLT block i (BLT_i), according to the dummy derivative
Expand All @@ -231,13 +239,16 @@ vConstraints[j] are the unknowns of eConstraints[j].
- `Brev::Vector{Int}`: Reverted equation association: `Erev[i] = der(E[k]) == E[i] ? k : 0`.
- `var_is_state::Function`: var_is_state(`vi`) returns true, if `vi` is defined to be a definite state.
- showMessage: showMessage2 below to print error/warning message
"""
function getConstraintSets(BLT_i::Vector{Int},
assignRev::Vector{Int},
Arev::Vector{Int},
Brev::Vector{Int},
var_is_state::Function,
showMessage::Function)
# Determine unknowns of BLT block i
BLT_i_unknowns = fill(0, length(BLT_i))
Expand Down Expand Up @@ -271,13 +282,25 @@ function getConstraintSets(BLT_i::Vector{Int},
pushfirst!(eConstraints, ceq) # move ceq to the beginning of the constraints vector

# Determine unknowns of ceq
veq = fill(0, 0)
veq = fill(0,0)
veq_states = fill(0,0)
for vc in vConstraints[1]
if Arev[vc] > 0
push!(veq, Arev[vc])
if var_is_state( Arev[vc] )
push!(veq_states, Arev[vc])
else
push!(veq, Arev[vc])
end
end
end
@assert( length(veq) > 0 )
if length(veq) < length(ceq)
nRemoveStates = length(veq_states) - (length(ceq) - length(veq))
showMessage("$nRemoveStates of the followig defined states cannot be states!";
severity = ERROR,
variables = veq_states,
equations = ceq)
return (eConstraints, vConstraints, false)
end
pushfirst!(vConstraints, veq) # move veq to the beginning of the constraints vector
end

Expand Down Expand Up @@ -479,7 +502,12 @@ mutable struct EquationGraph

# Initialize the active variables
vActive = fill(true, length(A))

for v = 1:length(vActive)
if Arev[v] > 0
vActive[ Arev[v] ] = false # Arev[v] is a potential state
end
end

# Define empty constraint sets
eConstraintsVec = Vector{Vector{Int}}[]
vConstraintsVec = Vector{Vector{Int}}[]
Expand Down Expand Up @@ -521,20 +549,21 @@ mutable struct EquationGraph
if highest
# Get all equation sets eConstraints and their corresponding unknowns vConstraints
# from lowest to highest differentiation order (eConstraints[end] is c)
(eConstraints, vConstraints, success2) = getConstraintSets(BLT_i, assignRev, Arev, Brev, showMessage)
(eConstraints, vConstraints, success2) = getConstraintSets(BLT_i, assignRev, Arev, Brev,
stateSelectionFunctions.var_is_state, showMessage)
if !success2
success=false
break
end

# Initialize vActive
for vc in vConstraints
for v in vc
if Arev[v] > 0
vActive[ Arev[v] ] = false # Arev[v] is a potential state
end
end
end
#for vc in vConstraints
# for v in vc
# if Arev[v] > 0
# vActive[ Arev[v] ] = false # Arev[v] is a potential state
# end
# end
#end

push!(eConstraintsVec, eConstraints)
push!(vConstraintsVec, vConstraints)
Expand Down Expand Up @@ -963,7 +992,7 @@ function sortEquations!(eq::EquationGraph)::Nothing
# Matching of all equations
assign = matching(eq.G, length(eq.A), eq.vActive)
assignRev = revertAssociation(assign)

# Sort equations and find strong components
blt = BLT(eq.G,assign)

Expand Down Expand Up @@ -1160,6 +1189,7 @@ function getSortedAndSolvedAST(G, # Typically ::Vector{Vector{Int}}
stateSelectionFunctions::StateSelectionFunctions;
log::Bool = false,
logDetails::Bool = false,
logStates::Bool = false,
modelName = "???",
unitless::Bool = false,
defaultParameterAndStartValues = nothing)
Expand Down Expand Up @@ -1270,16 +1300,29 @@ function getSortedAndSolvedAST(G, # Typically ::Vector{Vector{Int}}

# Select dummy states
if i < length(eConstraints)
# Not on highest derivative level. Solved variables are dummy states.
# Not on highest derivative level. Solved variables are dummy states
for v in eq.vSolved
@assert(!eq.vActive[v])
eq.vActive[v] = true
end
if log
println(" All other unknowns are solved and are dummy states.")
end

if length(eConstraints[i]) == length(vConstraints[i])
# N equations in N unknowns - tearing variables are dummy states
for v in eq.vTear
@assert(!eq.vActive[v])
eq.vActive[v] = true
end
if log
println(" All unknowns are dummy states.")
end

else
if log
println(" All solved unknowns are dummy states.")
end
end
elseif log
println(" All other unknowns are solved.")
println(" All unknowns are solved.")
end

# Further actions depending on type of equations
Expand Down Expand Up @@ -1379,7 +1422,7 @@ function getSortedAndSolvedAST(G, # Typically ::Vector{Vector{Int}}
v_stateCategory, v_length, v_unit, v_fixed2,
v_nominal, v_unbounded))
end
if log
if log || logStates
println("\nSelected ODE states: ")
x_table = ModiaBase.get_x_table(x_info)
show(stdout, x_table; allrows=true, allcols=true, summary=false, eltypes=false)
Expand Down

0 comments on commit 09a1c47

Please sign in to comment.