Unclear number of States in simple example

johhell · January 24, 2023, 7:45am

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.

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.

********************
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

baggepinnen · January 24, 2023, 7:49am

The confusion comes from ModelingToolkit confusing “variables” with “state” and “state variables”. There are three equations

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

zdenek_hurak · January 24, 2023, 9:01am

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.

johhell · January 24, 2023, 9:07am

Thanks. Clear so far.
I think there is no way to eliminate this additional linear equation.

help?> structural_simplify

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)

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

baggepinnen · January 24, 2023, 9:13am

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.

baggepinnen · January 24, 2023, 10:02am

The problem is that the system has a singular mass matrix E, i.e., it’s a descriptor system on the form

\begin{aligned} Eẋ &= 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

baggepinnen · January 24, 2023, 1:01pm

If you linearize between two points and find a minimal realization, you get what I think you are looking for:

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> 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 is available by calling named_ss instead of linearize

using ControlSystemsMTK
nsys = named_ss(rlc_model2, [source.V.u], [capacitor.p.v], allow_input_derivatives=true)
minreal(nsys)

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)