# Unclear number of States in simple example
**URL:** https://discourse.julialang.org/t/unclear-number-of-states-in-simple-example/93446
**Category:** New to Julia
**Tags:** modelingtoolkit
**Created:** [January 24, 2023, 7:45am UTC](https://discourse.julialang.org/t/unclear-number-of-states-in-simple-example/93446 "2023-01-24T07:45:26Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@johhell](https://discourse.julialang.org/u/johhell)
#### Post date: [January 24, 2023, 7:45am UTC](https://discourse.julialang.org/t/unclear-number-of-states-in-simple-example/93446/1 "2023-01-24T07:45:26Z")
I tested a simple component-based electrical example with a resistor, capacitor and an inductor. An additional resistor is connected in parallel with the inductor.
```julia
using ModelingToolkit
using LinearAlgebra
@variables t
D = Differential(t)
@connector function Pin(;name)
sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
ODESystem(Equation[], t, sts, []; name=name)
end
function Ground(;name)
@named g = Pin()
eqs = [g.v ~ 0]
compose(ODESystem(eqs, t, [], []; name=name), g)
end
function OnePort(;name)
@named p = Pin()
@named n = Pin()
sts = @variables v(t)=1.0 i(t)=1.0
eqs = [
v ~ p.v - n.v
0 ~ p.i + n.i
i ~ p.i ]
compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end
function Resistor(;name, R = 1.0)
@named oneport = OnePort()
@unpack v, i = oneport
ps = @parameters R=R
eqs = [v ~ i * R]
extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end
function Capacitor(;name, C = 1.0)
@named oneport = OnePort()
@unpack v, i = oneport
ps = @parameters C=C
eqs = [D(v) ~ i / C]
extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end
function ConstantVoltage(;name, V = 1.0)
@named oneport = OnePort()
@unpack v = oneport
ps = @parameters V=V
eqs = [V ~ v]
extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end
function Inductor(; name, L=1.0)
@named oneport = OnePort()
@unpack v, i = oneport
pars = @parameters L = L
eqs = [v ~ L * D(i)]
extend(ODESystem(eqs, t, [], pars; name = name), oneport)
end
R = 0.5; C = 0.2; L= 0.6; V = 1.0
@named resistor = Resistor(R=R)
@named resistor2 = Resistor(R=5.0)
@named capacitor = Capacitor(C=C)
@named inductor = Inductor(L=L)
@named source = ConstantVoltage(V=V)
@named ground = Ground()
connections2 = [
connect(source.p, resistor.p)
connect(resistor.n, inductor.p)
connect(inductor.n, capacitor.p)
connect(capacitor.n, source.n)
connect(capacitor.n, ground.g)
connect(inductor.p, resistor2.p) # resistor2 // inductor
connect(inductor.n, resistor2.n) # resistor2 // inductor
]
@named rlc_model2 = ODESystem(connections2, t, systems=[resistor, capacitor, source, ground, inductor, resistor2])
sys2 = structural_simplify(rlc_model2)
printstyled("\n ********************\nR - C - L//R2\n******************** \n", color=:red, bold=true)
@show states(sys2)
lsys2, ssys = linearize(rlc_model2, [], [])
printstyled("Linearization\n", underline=true)
show(IOContext(stdout, :compact => false), "text/plain", lsys2.A)
println()
printstyled("Eigenvalues\n", underline=true)
eigvals(lsys2.A)
```
I expected 2 states, but the result was different.
```julia
********************
R - C - L//R2
********************
states(sys2) = Any[capacitor₊v(t), inductor₊i(t), resistor2₊i(t)]
Linearization
3×3 Array{Float64, 2}:
0.0 5.0 5.0
-1.6666666666666667 -0.8333333333333334 -0.8333333333333334
0.15151515151515152 -0.8333333333333334 -0.8333333333333334
Eigenvalues
3-element Vector{ComplexF64}:
-0.8333333333333337 - 2.623225711087996im
-0.8333333333333337 + 2.623225711087996im
2.2624174687408427e-16 + 0.0im
```
remark: Test with `MODIA` gave the correct result.
Thanks / Johann
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [January 24, 2023, 7:49am UTC](https://discourse.julialang.org/t/unclear-number-of-states-in-simple-example/93446/2 "2023-01-24T07:49:31Z")
The confusion comes from ModelingToolkit confusing “variables” with “state” and “state variables”. There are three equations
```julia
julia> equations(sys2)
3-element Vector{Equation}:
Differential(t)(capacitor₊v(t)) ~ capacitor₊i(t) / capacitor₊C
Differential(t)(inductor₊i(t)) ~ resistor2₊v(t) / inductor₊L
0 ~ resistor2₊R*resistor2₊i(t) - resistor2₊v(t)
```
of which one is an algebraic equation. The system thus has two state variables just like you expected, and one algebraic variable. It’s unfortunate and confusing to a lot of people that ModelingToolkit, in contrast to the literature, call all of them states
---
### Author:  [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)
#### Post date: [January 24, 2023, 9:01am UTC](https://discourse.julialang.org/t/unclear-number-of-states-in-simple-example/93446/3 "2023-01-24T09:01:10Z")
> [@baggepinnen](#):
>
> The system thus has two states just like you expected, and one algebraic variable. It’s unfortunate and confusing to a lot of people that ModelingToolkit, in contrast to the literature, call all of them states
I do not want to be nitpicking, but here I would argue that even calling voltages and currents as **states** is unfortunate. They are **state variables**. It is their **values** (in some literature, especially on discrete-event systems and hybrid systems they even call it valuations) that determine the state at which the system is found.
The distinction is actually quite important because without it you can see/hear people trapped in trying to interpret concepts such as controllability and observability with respect to the individual (scalar) variables while these are defined with respect to the states (that is, vectors in a vector state space).
Similarly, saying that the system has just two states – the current through an inductor and a voltage on the capacitor – is confusing because the system can be actually found in an infinite number of states. The state space is a vector space of dimension two, but the number of elements of this states space – the number of states – is infinite.
---
### Author:  [@johhell](https://discourse.julialang.org/u/johhell)
#### Post date: [January 24, 2023, 9:07am UTC](https://discourse.julialang.org/t/unclear-number-of-states-in-simple-example/93446/4 "2023-01-24T09:07:28Z")
Thanks. Clear so far.
I think there is no way to eliminate this additional linear equation.
### help?\> structural\_simplify
```julia-auto
search: structural_simplify
structural_simplify(sys)
structural_simplify(sys, io; simplify, kwargs...)
Structurally simplify algebraic equations in a system and compute the topological sort of the
observed equations. When simplify=true, the simplify function will be applied during the
tearing process. It also takes kwargs allow_symbolic=false and allow_parameter=true which
limits the coefficient types during tearing.
The optional argument io may take a tuple (inputs, outputs). This will convert all inputs to
parameters and allow them to be unconnected, i.e., simplification will allow models where
n_states = n_equations - n_inputs.
```
### small modification
I tried to modify the example with `resistor2` in the first position of `systems`
`@named rlc_model2 = ODESystem(connections2, t, systems=[resistor2, resistor, capacitor, source, ground, inductor])`
The result of matrix `A` is different now, but EW’s are correct (with an additional unwanted value)
```julia-auto
states(sys2) = Any[capacitor₊v(t), inductor₊i(t), inductor₊p₊v(t)]
Linearization
3×3 Array{Float64, 2}:
-0.0 -0.0 -10.0
-0.0 -8.333333333333334 -16.666666666666668
0.0 3.7878787878787885 6.666666666666668
Eigenvalues
3-element Vector{ComplexF64}:
-0.8333333333333331 - 2.623225711087996im
-0.8333333333333331 + 2.623225711087996im
-0.0 + 0.0im
```
BR Johann
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [January 24, 2023, 9:13am UTC](https://discourse.julialang.org/t/unclear-number-of-states-in-simple-example/93446/5 "2023-01-24T09:13:17Z")
I agree, the wording I used was sloppy as well, the state would be composed out of state variables. I’ve updated the post to use a more accurate language.
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [January 24, 2023, 10:02am UTC](https://discourse.julialang.org/t/unclear-number-of-states-in-simple-example/93446/6 "2023-01-24T10:02:34Z")
> [@johhell](#):
>
> The result of matrix `A` is different now, but EW’s are correct (with an additional unwanted value)
The problem is that the system has a singular mass matrix E, i.e., it’s a descriptor system on the form
\begin{aligned} Eẋ &= Ax + Bu\\ y &= Cx + Du \end{aligned}
Currently, `linearize` does not output the mass matrix E but it probably should, I’ll submit a PR.
Once you have E, you can perform additional simplification using linear algebra. Ideally, though, MTK should be able to perform this simplification for linear cases like this.
* * *
Edit: Actually, `linearize` already performs the required simplification into a standard statespace system with non-singular mass matrix, but the resulting system requires input derivative information
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [January 24, 2023, 1:01pm UTC](https://discourse.julialang.org/t/unclear-number-of-states-in-simple-example/93446/7 "2023-01-24T13:01:12Z")
If you linearize between two points and find a minimal realization, you get what I think you are looking for:
```julia
using ModelingToolkit
using LinearAlgebra
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks
@variables t
D = Differential(t)
R = 0.5; C = 0.2; L= 0.6; V = 1.0
@named resistor = Resistor(R=R)
@named resistor2 = Resistor(R=5.0)
@named capacitor = Capacitor(C=C)
@named inductor = Inductor(L=L)
@named source = Voltage()
@named ground = Ground()
@named ref = Blocks.Constant(k=1)
connections2 = [
# connect(ref.output, :u, source.V) # Add this to simulate
connect(source.p, resistor.p)
connect(resistor.n, inductor.p)
connect(inductor.n, capacitor.p)
connect(capacitor.n, :y, source.n)
connect(capacitor.n, ground.g)
connect(inductor.p, resistor2.p)
connect(inductor.n, resistor2.n)
]
@named rlc_model2 = ODESystem(connections2, t, systems=[resistor, capacitor, source, ground, inductor, resistor2, ref])
lsys2, ssys = linearize(rlc_model2, [source.V.u], [capacitor.p.v], allow_input_derivatives=true)
using ControlSystemsBase
G = minreal(ss(lsys2...))
poles(G)
```
```julia
julia> G = minreal(ss(lsys2...))
ControlSystemsBase.StateSpace{Continuous, Float64}
A =
-1.666666666666666 -1.0713739108887084
7.0710678118654755 0.0
B =
-1.0713739108887081 -0.12856486930664499
0.0 0.0
C =
0.0 -1.0
D =
0.0 0.0
Continuous-time state-space model
julia> poles(G)
2-element Vector{ComplexF64}:
-0.833333333333333 + 2.623225711087998im
-0.833333333333333 - 2.623225711087998im
```
Note that this produces a model that requires input derivatives.
[Further simplification](https://juliacontrol.github.io/ControlSystemsMTK.jl/dev/#Internals:-Transformation-of-non-proper-models-to-proper-statespace-form) is available by calling [`named_ss`](https://juliacontrol.github.io/ControlSystemsMTK.jl/dev/api/#RobustAndOptimalControl.named_ss) instead of `linearize`
```julia
using ControlSystemsMTK
nsys = named_ss(rlc_model2, [source.V.u], [capacitor.p.v], allow_input_derivatives=true)
minreal(nsys)
```
```julia
julia> named_ss(rlc_model2, [source.V.u], [capacitor.p.v], allow_input_derivatives=true) |> minreal
NamedControlSystemsBase.StateSpace{Continuous, Float64}
A =
0.6039807824296497 2.4543162794335673
-3.6454898146825556 -2.2706474490963173
B =
-1.7919317184944736
0.0
C =
-0.5073245256547496 0.8437135739692423
D =
0.0
Continuous-time state-space model
With state names: x1 x2
input names: source₊V₊u(t)
output names: capacitor₊p₊v(t)
```