ModelingToolkit : Incorparating 2nd Order ODE's into functions

Specific Domains Modelling & Simulations
1

Hi all,

I am looking to incorporate the following valve functions in MTK from a paper the images below give the mathematical description of the valves. The paper can be found here

image
image
image

The way I am trying to implement the code is as follows

function OnePort(;name)
    @named in = Pin()
    @named out = Pin()
    sts = @variables Δp(t) = 0.0 q(t) = 0.0
    eqs = [
           Δp ~ out.p - in.p
           0 ~ in.q + out.q
           q ~ in.q
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), in, out)
end

function valvediff(Δp, AR, CQ)
    q = (Δp ≤ 0.0) *  CQ*AR*sqrt(sign(Δp)*Δp) +
        (Δp > 0.0) *  CQ*AR*sqrt(sign(Δp)*Δp)
return q
end


function ValveDiff(;name,CQ,Kp,Kf)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters CQ = CQ Kp = Kp Kf = Kf
    sts = @variables θ1(t) = 0.0 θ2(t) = 0.0 AR(t) = 0.0
    D = Differential(t)
    eqs = [
           D(θ1) ~ θ2    
           D(θ2) ~ (Δp ≤ 0.0) * Kp*Δp*cosd(θ1) - Kf*θ2  +
                   (Δp > 0.0) * Kp*Δp*cosd(θ1) - Kf*θ2

           AR ~ ((1-cosd(θ1))^2)/((1-cosd(90.0))^2)
           q ~ valvediff(Δp, AR, CQ)
          ]
    extend(ODESystem(eqs, t, sts, ps; name=name), oneport)
end

I then call this function by

@named AV =ValveDiff(CQ=CQ_AV, Kp = Kp_av, Kf = Kf_av)
@named MV =ValveDiff(CQ=CQ_MV, Kp = Kp_mv, Kf = Kf_mv)
@named TV =ValveDiff(CQ=CQ_TV, Kp = Kp_tv, Kf = Kf_tv)
@named PV =ValveDiff(CQ=CQ_PV, Kp = Kp_pv, Kf = Kf_pv)

Which is then used in the connect statement as follows

circ_eqs = [
    connect(LV.out, AV.in)
    connect(AV.out, SAS.in)
    connect(SAS.out, SAT.in)
    connect(SAT.out, SAR.in)
    connect(SAR.out, SCP.in)
    connect(SCP.out, SVN.in)
    connect(SVN.out, RA.in)
    connect(RA.out, TV.in)
    connect(TV.out, RV.in)
    connect(RV.out, PV.in)
    connect(PV.out, PAS.in)
    connect(PAS.out, PAT.in)
    connect(PAT.out, PAR.in)
    connect(PAR.out, PCP.in)
    connect(PCP.out, PVN.in)
    connect(PVN.out, LA.in)
    connect(LA.out, MV.in)
    connect(MV.out, LV.in)
]

The equations for the valves when when generated look as follows

 Differential(t)(AV₊θ1(t)) ~ AV₊θ2(t)
 Differential(t)(AV₊θ2(t)) ~ AV₊Kp*AV₊Δp(t)*cosd(AV₊θ1(t))*(AV₊Δp(t) <= 0.0) + AV₊Kp*AV₊Δp(t)*cosd(AV₊θ1(t))*(AV₊Δp(t) > 0.0) - (2//1)*AV₊Kf*AV₊θ2(t)
 Differential(t)(MV₊θ1(t)) ~ MV₊θ2(t)
 Differential(t)(MV₊θ2(t)) ~ MV₊Kp*MV₊Δp(t)*cosd(MV₊θ1(t))*(MV₊Δp(t) <= 0.0) + MV₊Kp*MV₊Δp(t)*cosd(MV₊θ1(t))*(MV₊Δp(t) > 0.0) - (2//1)*MV₊Kf*MV₊θ2(t)
 Differential(t)(PV₊θ1(t)) ~ PV₊θ2(t)
 Differential(t)(PV₊θ2(t)) ~ PV₊Kp*PV₊Δp(t)*cosd(PV₊θ1(t))*(PV₊Δp(t) <= 0.0) + PV₊Kp*PV₊Δp(t)*cosd(PV₊θ1(t))*(PV₊Δp(t) > 0.0) - (2//1)*PV₊Kf*PV₊θ2(t)
 Differential(t)(TV₊θ1(t)) ~ TV₊θ2(t)
 Differential(t)(TV₊θ2(t)) ~ TV₊Kp*TV₊Δp(t)*cosd(TV₊θ1(t))*(TV₊Δp(t) <= 0.0) + TV₊Kp*TV₊Δp(t)*cosd(TV₊θ1(t))*(TV₊Δp(t) > 0.0) - (2//1)*TV₊Kf*TV₊θ2(t)

When this example is done with a valve function I know which works, the solution correctly matches the literature so I am certain that the problem lies within the definition of the valve function somewhere.

I could set this up in the traditional way with a load of coupled ODEs however would like to use MTK for this problem. I have racked my brains for a couple of weeks now but I can’t think of any intuitive way of including the ODE into the function.

The above ODEs have transformed into first order ODEs.

Any help appreciated, if any other info is needed please reach out :slight_smile: - I posted the long message with detail as I think would serve well for others

Cheers,
Harry :slight_smile:

2

Those equations don’t look obviously wrong to me, what problem do you see when you run this? A fully runable (but perhaps simplified with fewer parts?) example that demonstrated the problem would help.

That said, these two expressions look redundant:

function valvediff(Δp, AR, CQ)
    q = (Δp ≤ 0.0) *  CQ*AR*sqrt(sign(Δp)*Δp) +
        (Δp > 0.0) *  CQ*AR*sqrt(sign(Δp)*Δp)
return q
end

Since the same quantity is evaluated independent of the sign of Δp, this could be just

valvediff(Δp, AR, CQ) = CQ*AR*sqrt(abs(Δp))

Similarly for

           D(θ2) ~ (Δp ≤ 0.0) * Kp*Δp*cosd(θ1) - Kf*θ2  +
                   (Δp > 0.0) * Kp*Δp*cosd(θ1) - Kf*θ2

The pressure sign doesn’t change anything, so that could be just

           D(θ2) ~ Kp*Δp*cosd(θ1) - Kf*θ2
3

The main problem is just extreme values and the solution diverging and I am using literature values for initial conditions however I could be missing something somewhere so i will have to check them in depth again.

As a heads up the IC I use for the valves are taken from the figure below

image

I will drop a MWE below with both a physiological valve function and the one above which is returning a non-physiological result.

using ModelingToolkit, DifferentialEquations, Plots 

### Parameter Values ###
τ = 0.85
Eₘᵢₙ = 0.03
Eₘₐₓ = 1.5
τes_lv = 0.3
τed_lv = 0.45

CQ_AV = 350.0
CQ_MV = 400.0
Kp_av = 5500.0
Kf_av = 50.0
Kf_mv = 50.0
Kp_mv = 5500.0

Rs = 1.11
Csa = 1.13
Csv = 11.0
MCFP = 7.0

### Start Modelling ###
@parameters t

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

function OnePort(;name)
    @named in = Pin()
    @named out = Pin()
    sts = @variables Δp(t) = 0.0 q(t) = 0.0
    eqs = [
           Δp ~ out.p - in.p
           0 ~ in.q + out.q
           q ~ in.q
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), in, out)
end

function Resistor(;name, R=1.0)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters R = R
    eqs = [
           Δp ~ - q * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C=1.0)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters C = C
    D = Differential(t)
    eqs = [
           D(Δp) ~ - q / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Ground(;name, P=0.0)
    @named g = Pin()
    ps = @parameters P = P
    eqs = [g.p ~ P]
    compose(ODESystem(eqs, t, [], ps; name=name), g)
end

function ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

    tᵢ = rem(t + (1 - Eshift) * τ, τ)

    Eₚ    = (tᵢ <= τₑₛ)               * ( 1 - cos(  tᵢ / τₑₛ         *pi))/2 +
            (tᵢ >  τₑₛ) * (tᵢ <= τₑₚ) * ( 1 + cos( (tᵢ-τₑₛ)/(τₑₚ-τₑₛ)*pi))/2 +
            (tᵢ <= τₑₚ)               * 0

    E     = Eₘᵢₙ + ( Eₘₐₓ -  Eₘᵢₙ ) * Eₚ

    return E
end

function DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

    tᵢ = rem(t + (1 - Eshift) * τ, τ)

    DEₚ  = (tᵢ <= τₑₛ)     * pi/ τₑₛ          * sin(tᵢ / τₑₛ         *pi)/2 +
            (tᵢ >  τₑₛ) * (tᵢ <= τₑₚ)  * pi/(τₑₚ - τₑₛ) * -sin((tᵢ-τₑₛ)/(τₑₚ - τₑₛ)         *pi)/2
            (tᵢ <= τₑₚ)               * 0
    DE   = ( Eₘₐₓ -  Eₘᵢₙ ) * DEₚ

    return DE
end

function ShiChamber(;name, V₀, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift=0.0, Ev=Inf)
    @named in = Pin()
    @named out = Pin()
    sts = @variables V(t) = 2.0 p(t) = 0.0
    ps = @parameters V₀ = V₀ Eₘᵢₙ = Eₘᵢₙ Eₘₐₓ = Eₘₐₓ τ = τ τₑₛ = τₑₛ τₑₚ = τₑₚ Eshift=Eshift Ev=Ev

    D = Differential(t)
    E  = ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)
    DE = DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

    eqs = [
        0 ~ in.p - out.p
        p ~ in.p
        V ~ p / E + V₀
        D(p) ~ (in.q + out.q) * E / (1 + 1/Ev * E) + p / (E * (1 + 1/Ev * E)) * DE
          ]

    compose(ODESystem(eqs, t, sts, ps; name=name), in, out)
end

#### Valve functions ###
function Valvesqrt(;name, CQ=1.0)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters CQ = CQ
    eqs = [
            q ~  (Δp < 0)*CQ*sqrt(sign(Δp)*Δp)
            ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end
## ODE valves
valvediff(Δp, AR, CQ) = CQ*AR*sqrt(abs(Δp))

function ValveDiff(;name,CQ,Kp,Kf)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters CQ = CQ Kp = Kp Kf = Kf
    sts = @variables θ1(t) = 0.0 θ2(t) = 0.0 AR(t) = 0.0
    D = Differential(t)
    θ_max = 75.0
    eqs = [
           D(θ1) ~ θ2    
           D(θ2) ~  Kp*sign(Δp)*Δp*cosd(θ1) - Kf*θ2  
           AR ~ ((1-cosd(θ1))^2)/((1-cosd(θ_max))^2)
           q ~ valvediff(Δp, AR, CQ)
          ]
    extend(ODESystem(eqs, t, sts, ps; name=name), oneport)
end

## Named Elements
@named LV = ShiChamber(V₀=0.0, τ=τ, Eₘₐₓ=Eₘₐₓ, Eₘᵢₙ=Eₘᵢₙ, τₑₛ=τes_lv, τₑₚ=τed_lv, Eshift=0.0, Ev=Inf)

## Changing to Valvesqrt is the working model 
@named AV = Valvesqrt(CQ=CQ_AV)#ValveDiff(CQ=CQ_AV, Kp = Kp_av, Kf = Kf_av)
@named MV = Valvesqrt(CQ=CQ_MV)#ValveDiff(CQ=CQ_MV, Kp = Kp_mv, Kf = Kf_mv)

@named Rs = Resistor(R=Rs)
@named Csa = Capacitor(C=Csa)
@named Csv = Capacitor(C=Csv)

@named ground = Ground() 

## Connect the system 
circ_eqs = [
    connect(LV.out, AV.in)
    connect(AV.out, Csa.in, Rs.in)
    connect(Rs.out, Csv.in, MV.in)
    connect(MV.out, LV.in)
    connect(Csa.out, ground.g)
    connect(Csv.out, ground.g)
    # Capacitor is connected to cavity
    # pressure which is same as 
    #C.out.p ~ 0.0
]


## Compose the whole ODE system
@named _circ_model = ODESystem(circ_eqs, t)

##
@named circ_model = compose(_circ_model,
                          [LV, AV, MV, Rs, Csa, Csv, ground])
                     


## And simplify it
circ_sys = structural_simplify(circ_model)

## Inital conditions
# Simple valve 
u0 = [MCFP, -MCFP, -MCFP]

# ODE valve 
#u0 = [MCFP, 75.0, 2.0, 0.0, 0.0, -MCFP, -MCFP]
prob = ODEProblem(circ_sys, u0, (0.0, 20.0))

@time sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-12)

#Plot
plot(sol, vars=[AV.q, MV.q], tspan=(17 * τ, 18 * τ))
plot(sol, vars=[LV.p, Csa.in.p, Csv.in.p], tspan=(16 * τ, 18 * τ), ylabel = "Pressure")
## Plot of how valve opening changes 
plot(sol, vars=[AV.θ1, MV.θ1], tspan=(17 * τ, 18 * τ))

You should be able to pick this up and run it immediately to give you a result for a working valve function, to change to the ODE valve we just need to change the AV and MV named components and then change the corresponding initial conditions.

I am not sure if it is how I implement the the second order or just a pesky negative somewhere :slight_smile:

Cheers,
Harry

4

The second-order stuff looks correct, but based on your hint about a sign error it is a bit suspicious that valves and resistors have opposite sign, is that intended? Your definition of Δp ~ out.p - in.p seem unusual to me, I would define positive pressure to produce positive flow (positive pressure implies in>out), but this definition has that reversed.

5

It is intentional I am building a whole library of things and some of things were defined in a negative way. Implemented the example below changing everything round and still getting extreme values as you will see

using ModelingToolkit, DifferentialEquations, Plots 

### Parameter Values ###
τ = 0.85
Eₘᵢₙ = 0.03
Eₘₐₓ = 1.5
τes_lv = 0.3
τed_lv = 0.45

CQ_AV = 350.0
CQ_MV = 400.0
Kp_av = 5500.0
Kf_av = 50.0
Kf_mv = 50.0
Kp_mv = 5500.0

Rs = 1.11
Csa = 1.13
Csv = 11.0
MCFP = 7.0

### Start Modelling ###
@parameters t

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

function OnePort(;name)
    @named in = Pin()
    @named out = Pin()
    sts = @variables Δp(t) = 0.0 q(t) = 0.0
    eqs = [
           Δp ~ in.p - out.p #out.p - in.p
           0 ~ in.q .+ out.q
           q ~ in.q
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), in, out)
end

function Resistor(;name, R=1.0)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters R = R
    eqs = [
           Δp ~  q * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C=1.0)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters C = C
    D = Differential(t)
    eqs = [
           D(Δp) ~  q / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Ground(;name, P=0.0)
    @named g = Pin()
    ps = @parameters P = P
    eqs = [g.p ~ P]
    compose(ODESystem(eqs, t, [], ps; name=name), g)
end

function ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

    tᵢ = rem(t + (1 - Eshift) * τ, τ)

    Eₚ    = (tᵢ <= τₑₛ)               * ( 1 - cos(  tᵢ / τₑₛ         *pi))/2 +
            (tᵢ >  τₑₛ) * (tᵢ <= τₑₚ) * ( 1 + cos( (tᵢ-τₑₛ)/(τₑₚ-τₑₛ)*pi))/2 +
            (tᵢ <= τₑₚ)               * 0

    E     = Eₘᵢₙ + ( Eₘₐₓ -  Eₘᵢₙ ) * Eₚ

    return E
end

function DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

    tᵢ = rem(t + (1 - Eshift) * τ, τ)

    DEₚ  = (tᵢ <= τₑₛ)     * pi/ τₑₛ          * sin(tᵢ / τₑₛ         *pi)/2 +
            (tᵢ >  τₑₛ) * (tᵢ <= τₑₚ)  * pi/(τₑₚ - τₑₛ) * -sin((tᵢ-τₑₛ)/(τₑₚ - τₑₛ)         *pi)/2
            (tᵢ <= τₑₚ)               * 0
    DE   = ( Eₘₐₓ -  Eₘᵢₙ ) * DEₚ

    return DE
end

function ShiChamber(;name, V₀, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift=0.0, Ev=Inf)
    @named in = Pin()
    @named out = Pin()
    sts = @variables V(t) = 2.0 p(t) = 0.0
    ps = @parameters V₀ = V₀ Eₘᵢₙ = Eₘᵢₙ Eₘₐₓ = Eₘₐₓ τ = τ τₑₛ = τₑₛ τₑₚ = τₑₚ Eshift=Eshift Ev=Ev

    D = Differential(t)
    E  = ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)
    DE = DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

    eqs = [
        0 ~ in.p - out.p
        p ~ in.p
        V ~ p / E + V₀
        D(p) ~ (in.q + out.q) * E / (1 + 1/Ev * E) + p / (E * (1 + 1/Ev * E)) * DE
          ]

    compose(ODESystem(eqs, t, sts, ps; name=name), in, out)
end

#### Valve functions ###
function Valvesqrt(;name, CQ=1.0)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters CQ = CQ
    eqs = [
            q ~  (Δp > 0)*CQ*sqrt(abs(Δp))
            ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

valvediff(Δp, AR, CQ) = CQ*AR*sqrt(abs(Δp))

function ValveDiff(;name,CQ,Kp,Kf)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters CQ = CQ Kp = Kp Kf = Kf
    sts = @variables θ1(t) = 0.0 θ2(t) = 0.0 AR(t) = 1.0
    D = Differential(t)
    θ_max = 75.0
    ## Comes from the Shi numerical simulation paper
    eqs = [
           D(θ1) ~ θ2    
           D(θ2) ~  Kp*Δp*cosd(θ1) - Kf*θ2  
           AR ~ ((1-cosd(θ1))^2)/((1-cosd(θ_max))^2)
           q ~ valvediff(Δp, AR, CQ)
          ]
    extend(ODESystem(eqs, t, sts, ps; name=name), oneport)
end
## Named Elements
@named LV = ShiChamber(V₀=0.0, τ=τ, Eₘₐₓ=Eₘₐₓ, Eₘᵢₙ=Eₘᵢₙ, τₑₛ=τes_lv, τₑₚ=τed_lv, Eshift=0.0, Ev=Inf)

## Changing to Valvesqrt is the working model 
@named AV = ValveDiff(CQ=CQ_AV, Kp = Kp_av, Kf = Kf_av)#Valvesqrt(CQ=CQ_AV)
@named MV = ValveDiff(CQ=CQ_MV, Kp = Kp_mv, Kf = Kf_mv)#Valvesqrt(CQ=CQ_MV)

@named Rs = Resistor(R=Rs)
@named Csa = Capacitor(C=Csa)
@named Csv = Capacitor(C=Csv)

@named ground = Ground() 

## Connect the system 
circ_eqs = [
    connect(LV.out, AV.in)
    connect(AV.out, Csa.in, Rs.in)
    connect(Rs.out, Csv.in, MV.in)
    connect(MV.out, LV.in)
    connect(Csa.out, ground.g)
    connect(Csv.out, ground.g)
    # Capacitor is connected to cavity
    # pressure which is same as 
    #C.out.p ~ 0.0
]


## Compose the whole ODE system
@named _circ_model = ODESystem(circ_eqs, t)

##
@named circ_model = compose(_circ_model,
                          [LV, AV, MV, Rs, Csa, Csv, ground])
                     


## And simplify it
circ_sys = structural_simplify(circ_model)

## Inital conditions
# Simple valve 
#u0 = [MCFP, MCFP, MCFP]

# ODE valve 
u0 = [MCFP, 74.0, 0.0, 0.0, 0.0, MCFP, MCFP]
prob = ODEProblem(circ_sys, u0, (0.0, 20.0))

@time sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-12)

#Plot
plot(sol, vars=[AV.q, MV.q], tspan=(17 * τ, 18 * τ))
plot(sol, vars=[LV.p, Csa.in.p, Csv.in.p], tspan=(6 * τ, 8 * τ), ylabel = "Pressure")
## Plot of how valve opening changes 
plot(sol, vars=[AV.θ1, MV.θ1], tspan=(17 * τ, 18 * τ))

This here is the solution vector element 2 and 4 represent the valve opening angle ideally I would like the angles to bounce between 0-75 but clearly that does not happen I would if I could be implementing it or if there is still a negative bouncing around the valve implementation hmmmm :confused:

image
image1205×240 16.9 KB

6

Changing the valve sign to match the resistor

valvediff(Δp, AR, CQ) = -sign(Δp)*CQ*AR*sqrt(abs(Δp))

Gives a nice plot
plot

7

Thanks for this however I quite embarrassingly can not reproduce what you have got above, could you drop the working code for me to compare perhaps?

All these years and I am still fooled by minus signs :face_with_diagonal_mouth:

Cheers :slight_smile:

8 
using ModelingToolkit, DifferentialEquations, Plots 

### Parameter Values ###
τ = 0.85
Eₘᵢₙ = 0.03
Eₘₐₓ = 1.5
τes_lv = 0.3
τed_lv = 0.45

CQ_AV = 350.0
CQ_MV = 400.0
Kp_av = 5500.0
Kf_av = 50.0
Kf_mv = 50.0
Kp_mv = 5500.0

Rs = 1.11
Csa = 1.13
Csv = 11.0
MCFP = 7.0

### Start Modelling ###
@parameters t

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

function OnePort(;name)
    @named in = Pin()
    @named out = Pin()
    sts = @variables Δp(t) = 0.0 q(t) = 0.0
    eqs = [
           Δp ~ out.p - in.p
           0 ~ in.q + out.q
           q ~ in.q
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), in, out)
end

function Resistor(;name, R=1.0)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters R = R
    eqs = [
           Δp ~ - q * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C=1.0)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters C = C
    D = Differential(t)
    eqs = [
           D(Δp) ~ - q / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Ground(;name, P=0.0)
    @named g = Pin()
    ps = @parameters P = P
    eqs = [g.p ~ P]
    compose(ODESystem(eqs, t, [], ps; name=name), g)
end

function ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

    tᵢ = rem(t + (1 - Eshift) * τ, τ)

    Eₚ    = (tᵢ <= τₑₛ)               * ( 1 - cos(  tᵢ / τₑₛ         *pi))/2 +
            (tᵢ >  τₑₛ) * (tᵢ <= τₑₚ) * ( 1 + cos( (tᵢ-τₑₛ)/(τₑₚ-τₑₛ)*pi))/2 +
            (tᵢ <= τₑₚ)               * 0

    E     = Eₘᵢₙ + ( Eₘₐₓ -  Eₘᵢₙ ) * Eₚ

    return E
end

function DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

    tᵢ = rem(t + (1 - Eshift) * τ, τ)

    DEₚ  = (tᵢ <= τₑₛ)     * pi/ τₑₛ          * sin(tᵢ / τₑₛ         *pi)/2 +
            (tᵢ >  τₑₛ) * (tᵢ <= τₑₚ)  * pi/(τₑₚ - τₑₛ) * -sin((tᵢ-τₑₛ)/(τₑₚ - τₑₛ)         *pi)/2
            (tᵢ <= τₑₚ)               * 0
    DE   = ( Eₘₐₓ -  Eₘᵢₙ ) * DEₚ

    return DE
end

function ShiChamber(;name, V₀, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift=0.0, Ev=Inf)
    @named in = Pin()
    @named out = Pin()
    sts = @variables V(t) = 2.0 p(t) = 0.0
    ps = @parameters V₀ = V₀ Eₘᵢₙ = Eₘᵢₙ Eₘₐₓ = Eₘₐₓ τ = τ τₑₛ = τₑₛ τₑₚ = τₑₚ Eshift=Eshift Ev=Ev

    D = Differential(t)
    E  = ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)
    DE = DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

    eqs = [
        0 ~ in.p - out.p
        p ~ in.p
        V ~ p / E + V₀
        D(p) ~ (in.q + out.q) * E / (1 + 1/Ev * E) + p / (E * (1 + 1/Ev * E)) * DE
          ]

    compose(ODESystem(eqs, t, sts, ps; name=name), in, out)
end

#### Valve functions ###
function Valvesqrt(;name, CQ=1.0)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters CQ = CQ
    eqs = [
            q ~  (Δp < 0)*CQ*sqrt(sign(Δp)*Δp)
            ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end
## ODE valves
# I'm pretty sure this is the only line I changed.
valvediff(Δp, AR, CQ) = -sign(Δp)*CQ*AR*sqrt(abs(Δp))

function ValveDiff(;name,CQ,Kp,Kf)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters CQ = CQ Kp = Kp Kf = Kf
    sts = @variables θ1(t) = 0.0 θ2(t) = 0.0 AR(t) = 0.0
    D = Differential(t)
    θ_max = 75.0
    eqs = [
           D(θ1) ~ θ2    
           D(θ2) ~  Kp*sign(Δp)*Δp*cosd(θ1) - Kf*θ2  
           AR ~ ((1-cosd(θ1))^2)/((1-cosd(θ_max))^2)
           q ~ valvediff(Δp, AR, CQ)
          ]
    extend(ODESystem(eqs, t, sts, ps; name=name), oneport)
end

## Named Elements
@named LV = ShiChamber(V₀=0.0, τ=τ, Eₘₐₓ=Eₘₐₓ, Eₘᵢₙ=Eₘᵢₙ, τₑₛ=τes_lv, τₑₚ=τed_lv, Eshift=0.0, Ev=Inf)

## Changing to Valvesqrt is the working model 
@named AV = Valvesqrt(CQ=CQ_AV)#ValveDiff(CQ=CQ_AV, Kp = Kp_av, Kf = Kf_av)
@named MV = Valvesqrt(CQ=CQ_MV)#ValveDiff(CQ=CQ_MV, Kp = Kp_mv, Kf = Kf_mv)

@named Rs = Resistor(R=Rs)
@named Csa = Capacitor(C=Csa)
@named Csv = Capacitor(C=Csv)

@named ground = Ground() 

## Connect the system 
circ_eqs = [
    connect(LV.out, AV.in)
    connect(AV.out, Csa.in, Rs.in)
    connect(Rs.out, Csv.in, MV.in)
    connect(MV.out, LV.in)
    connect(Csa.out, ground.g)
    connect(Csv.out, ground.g)
    # Capacitor is connected to cavity
    # pressure which is same as 
    #C.out.p ~ 0.0
]


## Compose the whole ODE system
@named _circ_model = ODESystem(circ_eqs, t)

##
@named circ_model = compose(_circ_model,
                          [LV, AV, MV, Rs, Csa, Csv, ground])



## And simplify it
circ_sys = structural_simplify(circ_model)

## Inital conditions
# Simple valve 
u0 = [MCFP, -MCFP, -MCFP]

# ODE valve 
#u0 = [MCFP, 75.0, 2.0, 0.0, 0.0, -MCFP, -MCFP]
prob = ODEProblem(circ_sys, u0, (0.0, 20.0))

@time sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-12)

#Plot
p1 = plot(sol, vars=[AV.q, MV.q], tspan=(17 * τ, 18 * τ))
p2 = plot(sol, vars=[LV.p, Csa.in.p, Csv.in.p], tspan=(16 * τ, 18 * τ), ylabel = "Pressure")
## Plot of how valve opening changes 
# p3 = plot(sol, vars=[AV.θ1, MV.θ1], tspan=(17 * τ, 18 * τ))
plot(p1, p2)
``
9

Thanks for this,

However I know this works the function Valvesqrt I wrote myself

It is when i run the function ValveDiff for the AV and MV where I run into issues I will drop the example below which when run will give non-physiological results

The main issues are surrounding the function Valvediff and valvediff

using ModelingToolkit, DifferentialEquations, Plots 

### Parameter Values ###
τ = 0.85
Eₘᵢₙ = 0.03
Eₘₐₓ = 1.5
τes_lv = 0.3
τed_lv = 0.45

CQ_AV = 350.0
CQ_MV = 400.0
Kp_av = 5500.0
Kf_av = 50.0
Kf_mv = 50.0
Kp_mv = 5500.0

Rs = 1.11
Csa = 1.13
Csv = 11.0
MCFP = 7.0

### Start Modelling ###
@parameters t

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

function OnePort(;name)
    @named in = Pin()
    @named out = Pin()
    sts = @variables Δp(t) = 0.0 q(t) = 0.0
    eqs = [
           Δp ~ out.p - in.p
           0 ~ in.q + out.q
           q ~ in.q
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), in, out)
end

function Resistor(;name, R=1.0)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters R = R
    eqs = [
           Δp ~ - q * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C=1.0)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters C = C
    D = Differential(t)
    eqs = [
           D(Δp) ~ - q / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Ground(;name, P=0.0)
    @named g = Pin()
    ps = @parameters P = P
    eqs = [g.p ~ P]
    compose(ODESystem(eqs, t, [], ps; name=name), g)
end

function ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

    tᵢ = rem(t + (1 - Eshift) * τ, τ)

    Eₚ    = (tᵢ <= τₑₛ)               * ( 1 - cos(  tᵢ / τₑₛ         *pi))/2 +
            (tᵢ >  τₑₛ) * (tᵢ <= τₑₚ) * ( 1 + cos( (tᵢ-τₑₛ)/(τₑₚ-τₑₛ)*pi))/2 +
            (tᵢ <= τₑₚ)               * 0

    E     = Eₘᵢₙ + ( Eₘₐₓ -  Eₘᵢₙ ) * Eₚ

    return E
end

function DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

    tᵢ = rem(t + (1 - Eshift) * τ, τ)

    DEₚ  = (tᵢ <= τₑₛ)     * pi/ τₑₛ          * sin(tᵢ / τₑₛ         *pi)/2 +
            (tᵢ >  τₑₛ) * (tᵢ <= τₑₚ)  * pi/(τₑₚ - τₑₛ) * -sin((tᵢ-τₑₛ)/(τₑₚ - τₑₛ)         *pi)/2
            (tᵢ <= τₑₚ)               * 0
    DE   = ( Eₘₐₓ -  Eₘᵢₙ ) * DEₚ

    return DE
end

function ShiChamber(;name, V₀, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift=0.0, Ev=Inf)
    @named in = Pin()
    @named out = Pin()
    sts = @variables V(t) = 2.0 p(t) = 0.0
    ps = @parameters V₀ = V₀ Eₘᵢₙ = Eₘᵢₙ Eₘₐₓ = Eₘₐₓ τ = τ τₑₛ = τₑₛ τₑₚ = τₑₚ Eshift=Eshift Ev=Ev

    D = Differential(t)
    E  = ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)
    DE = DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

    eqs = [
        0 ~ in.p - out.p
        p ~ in.p
        V ~ p / E + V₀
        D(p) ~ (in.q + out.q) * E / (1 + 1/Ev * E) + p / (E * (1 + 1/Ev * E)) * DE
          ]

    compose(ODESystem(eqs, t, sts, ps; name=name), in, out)
end

#### Valve functions ###
## ODE valves
# I'm pretty sure this is the only line I changed.
valvediff(Δp, AR, CQ) = -sign(Δp)*CQ*AR*sqrt(abs(Δp))

function ValveDiff(;name,CQ,Kp,Kf)
    @named oneport = OnePort()
    @unpack Δp, q = oneport
    ps = @parameters CQ = CQ Kp = Kp Kf = Kf
    sts = @variables θ1(t) = 0.0 θ2(t) = 0.0 AR(t) = 0.0
    D = Differential(t)
    θ_max = 75.0
    eqs = [
           D(θ1) ~ θ2    
           D(θ2) ~  Kp*Δp*cosd(θ1) - Kf*θ2
           AR ~ ((1-cosd(θ1))^2)/((1-cosd(θ_max))^2)
           q ~ valvediff(Δp, AR, CQ)
          ]
    extend(ODESystem(eqs, t, sts, ps; name=name), oneport)
end
((1-cosd(0))^2)/((1-cosd(90))^2)
## Named Elements
@named LV = ShiChamber(V₀=0.0, τ=τ, Eₘₐₓ=Eₘₐₓ, Eₘᵢₙ=Eₘᵢₙ, τₑₛ=τes_lv, τₑₚ=τed_lv, Eshift=0.0, Ev=Inf)

@named AV = ValveDiff(CQ=CQ_AV, Kp = Kp_av, Kf = Kf_av)
@named MV = ValveDiff(CQ=CQ_MV, Kp = Kp_mv, Kf = Kf_mv)

@named Rs = Resistor(R=Rs)
@named Csa = Capacitor(C=Csa)
@named Csv = Capacitor(C=Csv)

@named ground = Ground() 

## Connect the system 
circ_eqs = [
    connect(LV.out, AV.in)
    connect(AV.out, Csa.in, Rs.in)
    connect(Rs.out, Csv.in, MV.in)
    connect(MV.out, LV.in)
    connect(Csa.out, ground.g)
    connect(Csv.out, ground.g)
    # Capacitor is connected to cavity
    # pressure which is same as 
    #C.out.p ~ 0.0
]


## Compose the whole ODE system
@named _circ_model = ODESystem(circ_eqs, t)

##
@named circ_model = compose(_circ_model,
                          [LV, AV, MV, Rs, Csa, Csv, ground])



## And simplify it
circ_sys = structural_simplify(circ_model)

u0 = [MCFP, 75.0, 0.0, 0.0, 0.0, -MCFP, -MCFP]
prob = ODEProblem(circ_sys, u0, (0.0, 20.0))

@time sol = solve(prob, reltol=1e-12, abstol=1e-12)

#Plot
p1 = plot(sol, vars=[AV.q, MV.q], tspan=(6 * τ, 7 * τ))
p2 = plot(sol, vars=[LV.p, Csa.in.p, Csv.in.p], tspan=(6 * τ, 7 * τ), ylabel = "Pressure")
p3 = plot(sol, vars=[AV.θ1, MV.θ1], tspan=(17 * τ, 18 * τ))
plot(p1, p2)

The code above now only contains the non-working example where the above contained both the working and non working example. I am still pushing on my side so I think it has some worth sharing this. I feel like the error has to be narrowed down to a missing minus sign sitting in the function definition

Sorry if there was any confusion and thank you for your help :slight_smile:

Cheers,
Harry