Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Dropped unknowns with non-fully determined structural simplifications #2515

Closed
ChrisRackauckas opened this issue Mar 1, 2024 · 1 comment · Fixed by #2587
Closed

Dropped unknowns with non-fully determined structural simplifications #2515

ChrisRackauckas opened this issue Mar 1, 2024 · 1 comment · Fixed by #2587
Assignees

Comments

@ChrisRackauckas
Copy link
Member

ChrisRackauckas commented Mar 1, 2024

The case of 25 equations with 26 unknowns is an undeterdetermined system which is then simplified down to 0 equations which should then have 1 unknown, since there's 25 observed equations, but some shortcut is making it turn into 0 unknowns.

using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D

@connector function Pin(; name)
    sts = @variables v(t) [guess = 1.0] i(t) [guess = 1.0, connect = Flow]
    ODESystem(Equation[], t, sts, []; name = name)
end

@component function Ground(; name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name = name), g)
end

@component function OnePort(; name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t) [guess = 1.0] i(t) [guess = 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

@component 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

@component 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

@component 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

@component function Inductor(; name, L = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters L = L
    eqs = [
        D(i) ~ v / L
    ]
    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
end

R = 1.0
C = 1.0
V = 1.0
@named resistor = Resistor(R = R)
@named capacitor = Capacitor(C = C)
@named source = ConstantVoltage(V = V)
@named ground = Ground()

rc_eqs = [connect(source.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source.n)
          connect(capacitor.n, ground.g)]

@named rc_model = ODESystem(rc_eqs, t)
rc_model = compose(rc_model, [resistor, capacitor, source, ground])

@named resistor2 = Resistor(R = R)
rc_eqs2 = [connect(source.p, resistor.p)
           connect(resistor.n, resistor2.p)
           connect(resistor2.n, capacitor.p)
           connect(capacitor.n, source.n)
           connect(capacitor.n, ground.g)]

@named _rc_model2 = ODESystem(rc_eqs2, t)
@named rc_model2 = compose(_rc_model2,
    [resistor, resistor2, capacitor, source, ground])
sys2 = structural_simplify(rc_model2)

u0 = [capacitor.v => 0.0
      capacitor.p.i => 0.0
      resistor.v => 0.0]
      
prob2 = ODEProblem(sys2, [], (0, 10.0), guesses = u0)
sol2 = solve(prob2, Rosenbrock23()) # Error

isys = generate_initializesystem(sys2)
isys = complete(isys) # 25 equations, 26 unknowns
iprob = NonlinearLeastSquaresProblem(isys, u0, [t=>0.0])
iprob[sys2.resistor.n.v]

isys2 = structural_simplify(isys; fully_determined = false) # 0 equations, 0 unknowns?
iprob2 = NonlinearLeastSquaresProblem(isys2, u0, [t=>0.0])
iprob2[sys2.resistor.n.v] # Error
ChrisRackauckas added a commit that referenced this issue Mar 1, 2024
@ChrisRackauckas
Copy link
Member Author

https://github.com/SciML/ModelingToolkit.jl/blob/master/test/components.jl#L88-L104

u0 = [capacitor.v => 0.0
      capacitor.p.i => 0.0
      resistor.v => 0.0]

IIRC just remove the last two.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants