I am working on ModelingToolkit.jl using the RC circuit example at https://mtk.sciml.ai/stable/tutorials/acausal_components/. It appears that as long as the Differential(t) terms in the equations are on the Left Hand Side, things work as expected. However, when the left and right hand sides of an equation are swapped (with no other changes), there is some trouble with finding the right state variables after alias elimination. Details are provided below:
I started with the flattened system of 20 equations and 20 variables:
eqs = [0 ~ resistor₊p₊i + source₊p₊i,
source₊p₊v ~ resistor₊p₊v,
0 ~ capacitor₊p₊i + resistor₊n₊i,
resistor₊n₊v ~ capacitor₊p₊v,
0 ~ capacitor₊n₊i + ground₊g₊i + source₊n₊i,
capacitor₊n₊v ~ source₊n₊v,
source₊n₊v ~ ground₊g₊v,
resistor₊v ~ resistor₊p₊v - resistor₊n₊v,
0 ~ resistor₊n₊i + resistor₊p₊i,
resistor₊i ~ resistor₊p₊i,
resistor₊v ~ resistorR*resistor₊i,
capacitor₊v ~ capacitor₊p₊v - capacitor₊n₊v,
0 ~ capacitor₊n₊i + capacitor₊p₊i,
capacitor₊i ~ capacitor₊p₊i,
Differential(t)(capacitor₊v) ~ capacitor₊i*(capacitorC^-1),
source₊v ~ source₊p₊v - source₊n₊v,
0 ~ source₊n₊i + source₊p₊i,
source₊i ~ source₊p₊i,
sourceV ~ source₊v,
ground₊g₊v ~ 0]
I created an ODESystem:
@named rc_model = ODESystem(eqs,t)
I did alias elimination on the system:
sys_alias = alias_elimination(rc_model)
Differential(t)(capacitor₊v(t)) ~ resistor₊i(t) / capacitorC
0 ~ source₊v(t) - capacitor₊v(t) - resistor₊v(t)
0 ~ resistorR*resistor₊i(t) - resistor₊v(t)
0 ~ source₊v(t) - sourceV
Following are the state variables after alias elimination:
states(sys_alias)
capacitor₊v(t)
resistor₊v(t)
resistor₊i(t)
source₊v(t)
I then performed structural simplification on the system:
sys = structural_simplify(rc_model)
Following are now the equations:
Differential(t)(capacitor₊v(t)) ~ resistor₊i(t) / capacitorC
0 ~ resistorR*resistor₊i(t) + capacitor₊v(t) - sourceV
And the states are:
states(sys)
capacitor₊v(t)
resistor₊i(t)
So far everything as worked perfectly and as expected.
Now observe what happens when I swap the left and right sides of the differential equation, all other equations remaining exactly the same as before:
capacitor₊i*(capacitorC^-1) ~ Differential(t)(capacitor₊v)
Here is the full equation set:
eqs2 = [0 ~ resistor₊p₊i + source₊p₊i,
source₊p₊v ~ resistor₊p₊v,
0 ~ capacitor₊p₊i + resistor₊n₊i,
resistor₊n₊v ~ capacitor₊p₊v,
0 ~ capacitor₊n₊i + ground₊g₊i + source₊n₊i,
capacitor₊n₊v ~ source₊n₊v,
source₊n₊v ~ ground₊g₊v,
resistor₊v ~ resistor₊p₊v - resistor₊n₊v,
0 ~ resistor₊n₊i + resistor₊p₊i,
resistor₊i ~ resistor₊p₊i,
resistor₊v ~ resistorR*resistor₊i,
capacitor₊v ~ capacitor₊p₊v - capacitor₊n₊v,
0 ~ capacitor₊n₊i + capacitor₊p₊i,
capacitor₊i ~ capacitor₊p₊i,
capacitor₊i*(capacitorC^-1) ~ Differential(t)(capacitor₊v),
source₊v ~ source₊p₊v - source₊n₊v,
0 ~ source₊n₊i + source₊p₊i,
source₊i ~ source₊p₊i,
sourceV ~ source₊v,
ground₊g₊v ~ 0]
I created another model with this equation set:
@named rc_model2 = ODESystem(eqs2,t)
Now I perform alias elimination:
sys2_alias = alias_elimination(rc_model2)
Here are the resulting equations:
0 ~ source₊v(t) - capacitor₊v(t) - resistor₊v(t)
0 ~ resistorR*resistor₊i(t) - resistor₊v(t)
capacitor₊i(t) / capacitorC ~ Differential(t)(capacitor₊v(t))
0 ~ source₊v(t) - sourceV
And here are the states:
states(sys2_alias)
capacitor₊v(t)
resistor₊v(t)
resistor₊i(t)
source₊v(t)
Observe how there are only four state variables identified, however there is one more variable in the equations , capacitor₊i(t).
Now when I perform structural simplification, it throws an error:
structural_simplify(sys2_alias)
ExtraVariablesSystemException: The system is unbalanced. There are 5 highest order derivative variables and 4 equations.
More variables than equations, here are the potential extra variable(s):
capacitor₊v(t)
capacitor₊i(t)
Stacktrace:
[1] error_reporting(sys::ODESystem, bad_idxs::Vector{Int64}, n_highest_vars::Int64, iseqs::Bool)
@ ModelingToolkit.StructuralTransformations C:\Users\Sushil\.julia\packages\ModelingToolkit\eSpDL\src\structural_transformation\utils.jl:72
[2] check_consistency(sys::ODESystem)
@ ModelingToolkit.StructuralTransformations C:\Users\Sushil\.julia\packages\ModelingToolkit\eSpDL\src\structural_transformation\utils.jl:103
[3] structural_simplify(sys::ODESystem; simplify::Bool)
@ ModelingToolkit C:\Users\Sushil\.julia\packages\ModelingToolkit\eSpDL\src\systems\abstractsystem.jl:815
[4] structural_simplify(sys::ODESystem)
@ ModelingToolkit C:\Users\Sushil\.julia\packages\ModelingToolkit\eSpDL\src\systems\abstractsystem.jl:814
[5] top-level scope
@ In[82]:1
[6] eval
@ .\boot.jl:360 [inlined]
[7] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:1116
The variables resistor₊i(t) and capacitor₊i(t) are redundant. If I manually replace one of these variables with the other, then it works - there are 4 equations and 4 variables.
So the question is: why does alias_elimination not correctly identify the state variables when the Differential(t) term is on the right side? Am I missing something?
Thanks for any feedback and comments in resolving this puzzle. ModelingToolkit.jl is a great package and so far it has been of great help in my work.
Below is the full code which throws error:
using ModelingToolkit
@parameters t resistorR capacitorC sourceV
@variables resistor₊p₊i(t) source₊p₊i(t) source₊p₊v(t) resistor₊p₊v(t) capacitor₊p₊i(t) resistor₊n₊i(t)
@variables resistor₊n₊v(t) capacitor₊p₊v(t) capacitor₊n₊i(t) ground₊g₊i(t) source₊n₊i(t) source₊v(t) source₊i(t)
@variables capacitor₊n₊v(t) source₊n₊v(t) ground₊g₊v(t) resistor₊v(t) resistor₊i(t) capacitor₊v(t) capacitor₊i(t)
eqs2 = [0 ~ resistor₊p₊i + source₊p₊i,
source₊p₊v ~ resistor₊p₊v,
0 ~ capacitor₊p₊i + resistor₊n₊i,
resistor₊n₊v ~ capacitor₊p₊v,
0 ~ capacitor₊n₊i + ground₊g₊i + source₊n₊i,
capacitor₊n₊v ~ source₊n₊v,
source₊n₊v ~ ground₊g₊v,
resistor₊v ~ resistor₊p₊v - resistor₊n₊v,
0 ~ resistor₊n₊i + resistor₊p₊i,
resistor₊i ~ resistor₊p₊i,
resistor₊v ~ resistorR*resistor₊i,
capacitor₊v ~ capacitor₊p₊v - capacitor₊n₊v,
0 ~ capacitor₊n₊i + capacitor₊p₊i,
capacitor₊i ~ capacitor₊p₊i,
capacitor₊i*(capacitorC^-1) ~ Differential(t)(capacitor₊v),
source₊v ~ source₊p₊v - source₊n₊v,
0 ~ source₊n₊i + source₊p₊i,
source₊i ~ source₊p₊i,
sourceV ~ source₊v,
ground₊g₊v ~ 0]
@named rc_model2 = ODESystem(eqs2,t)
sys2_alias = alias_elimination(rc_model2)
structural_simplify(sys2_alias)
. Hopefully it’s done in the next month.