Skip to content

structural_simplify bug (I think) #1178

@ordicker

Description

@ordicker

Hi,
I'm translating Simulink model into MTK and it is great! On Simulink it's about 10 seconds for 10 seconds simulation and in MTK it's about 250ms(!).

However I think there is a bug in structural_simplify. The following code doesn't work

using ModelingToolkit, Plots, DifferentialEquations

function gimbal_plant(;name, J=1.0, k_s=1.0, k_v=1.0, ω_brk=1.0, T_c=1.0)
    @variables t T(t) ω_rel(t) θ_rel(t) ω_body(t) ω(t)
    @parameters _J _k_s _k_v _ω_brk _T_c
    D = Differential(t)
    eqs = [
        D(θ_rel) ~ ω_rel
        ω ~ ω_rel + ω_body
        D(ω) ~ (T-_k_s*θ_rel-_k_v*ω_rel-_T_c*tanh(10*ω_rel/_ω_brk))/_J
    ]
    ODESystem(eqs, t; name, defaults==>0.0, ω_body=>0.0, ω_rel=>0.0,
                                      θ_rel=>0.0, T=>0.0,
                                      _T_c=>T_c, _ω_brk=>ω_brk,
                                      _J=>J, _k_s=>k_s, _k_v=>k_v])
end

function PID_factory(;name, kp=1.0, ki=0.0, kd=0.0)
    @parameters _kp _ki _kd
    @variables t e(t) Ie(t) de(t) u(t)
    D = Differential(t)
    
    eqs = [
        D(Ie) ~ e
        D(e) ~ de
        u ~ _ki*Ie+_kp*e+_kd*de
    ]
    ODESystem(eqs, t ;name, defaults=[u=>0.0,
                                      e=>0.0, de=>0.0, Ie=>0.0,
                                      _kp=>kp, _ki=>ki, _kd=>kd])
end

function lead_lag_factory(;name,k=1.0, zero=1.0, pole=0.0)
    @variables t input(t) dinput(t) output(t) doutput(t)
    @parameters _zero _pole _k
    D = Differential(t)
    eqs = [
        D(input) ~ dinput
        D(output) ~ doutput
        _pole*doutput + output ~ _k*(_zero*dinput + input)
    ]
    ODESystem(eqs, t ;name, defaults=[output=>0.0, doutput=>0.0,
                                      input=>0.0, dinput=>0.0,
                                      _zero=>zero, _pole=>pole, _k=>k])
end

function sensor(;name, bw=1)
    @variables t input(t) output(t) doutput(t) ddoutput(t)
    @parameters a2 a1
    D = Differential(t)
    eqs = [
        D(output) ~ doutput
        D(doutput) ~ ddoutput
        input ~ output+a1*doutput+a2*ddoutput
        #output ~ input-(4e-4*doutput+1e-6*ddoutput)
    ]
    ODESystem(eqs, t ;name, defaults=[output=>0.0, doutput=>0.0,
                                      ddoutput=>0.0, input=>0.0,
                                      a1=>2*0.7/(2π*bw), a2=>1/(2π*bw^2)])
end


function gimbal_test()
    # plant params
    J = 0.029
    k_s = 0.2
    k_v = 0.545
    ω_brk = 0.074
    T_c = 0.106
    T_brk = 0.177
    # controller params
    K_dc = 6878.7
    Z1 = 0.01
    Z2 = 0.029
    P2 = 0.16
    Z3 = 0.06
    P3 = 0.009
    Z4 = 0.01
    bw = 465.0

    @variables t
    @named pid1 = PID_factory(kp=Z1,ki=1.0,kd=0.0)
    @named lag = lead_lag_factory(k=1.0, zero=Z2, pole=P2)
    @named lead = lead_lag_factory(k=1.0, zero=Z3, pole=P3)
    @named pid2 = PID_factory(kp=Z4,ki=1.0,kd=0.0)
    @named p = gimbal_plant(J=J,k_s=k_s,k_v=k_v,ω_brk=ω_brk,T_c=T_c)
    @named gyro = sensor(bw=bw)
    connections = [
        p.ω_body ~ 16*π^2/180*sin(2π*2t)
        gyro.input ~ p.ω
        pid1.e ~ -(180/π)*gyro.output
        lead.input ~ pid1.u
        lag.input ~ lead.output 
        pid2.e ~ lag.output
        p.T ~ 0.339*min(max(K_dc*pid2.u,-2.5),2.5)]    
    
    @named connected = compose(ODESystem(connections,t),
                               pid1, pid2, lead, lag, p, gyro)

    sys = structural_simplify(connected)
    #sys = alias_elimination(connected)
    prob = ODEProblem(sys,[],(0.0,10.0))
    return prob, p
end


prob, p = gimbal_test()

with this error:

ERROR: 

... Internal error in Tearing.jl: vs[1] = 0.

Stacktrace:

If I'm using alias_elimination instead of structural_simplify the code works and really fast(compared to Simulink)

Is it a bug? Am I using it in the wrong way?

Metadata

Metadata

Assignees

Labels

No labels
No labels

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions