ModelingToolkit.jl: Handling Output/Input Dependency in a DAE Model

AndhikaTama · February 22, 2025, 3:32am

Hello everyone,

I’m currently working on a dynamic of microgas turbine system model in Julia using a system of Differential-Algebraic Equations (DAE).

I want to adapt modelica acausal modelling problem into Julia using modellingtoolkit.jl and i already create a code for the system that look like this :

using ModelingToolkit, DifferentialEquations, CoolProp, Plots
@independent_variables t
D = Differential(t)
const fluid = "Air"
using ModelingToolkit, DifferentialEquations, CoolProp, Plots
@independent_variables t
D = Differential(t)
const fluid = "Air"
PropsSI(out::AbstractString, name1::AbstractString, value1::Real, name2::AbstractString, value2::Real, fluid::AbstractString) = CoolProp.PropsSI(out, name1, value1, name2, value2, fluid)
@register_symbolic PropsSI(out::AbstractString, name1::AbstractString, value1::Real, name2::AbstractString, value2::Real, fluid::AbstractString)
function h(p_in, T_in)
    h = PropsSI("H", "P", p_in, "T", T_in, fluid)
    return PropsSI("H", "P", p_in, "T", T_in, fluid)
end
function hf(p_in, T_in)
    hf = PropsSI("H", "P", p_in, "T", T_in, fluid)
    return PropsSI("H", "P", p_in, "T", T_in, fluid)
end
function ρ(p_in, T_in)
    ρ = PropsSI("D", "P", p_in, "T", T_in, fluid)
    return PropsSI("D", "P", p_in, "T", T_in, fluid)
end
function u(p_in, T_in)
    u = PropsSI("D", "P", p_in, "T", T_in, fluid)
    return PropsSI("D", "P", p_in, "T", T_in, fluid)
end
function cv(T)
	cv = PropsSI("Cvmass", "T", T, "P", 101325, fluid)
	return PropsSI("Cvmass", "T", T, "P", 101325, fluid)
end
function cp(T)
	cp = PropsSI("Cpmass", "T", T, "P", 101325, fluid)
	return PropsSI("Cpmass", "T", T, "P", 101325, fluid)
end
@register_symbolic h(p_in, T_in)
@register_symbolic hf(p_in, T_in)
@register_symbolic ρ(p_in, T_in)
@register_symbolic u(p_in, T_in)
@register_symbolic cv(T)
@register_symbolic cp(T)
@connector  function Port(;name) 
vars = @variables begin 
    p(t), [input = true]
	T(t), [input = true]
	P(t)
    mdot(t), [connect = Flow]
end
ODESystem(Equation[], t, vars, [];name=name)
end
function Source(;name,source_temperature, source_pressure,source_mdot, source_power) 
    @named port = Port()
    para = @parameters begin
        
    end
    vars = @variables begin

     end

    eqs = [
        port.mdot ~ source_mdot 
		port.P ~ source_power
        port.p ~ source_pressure
		port.T ~ source_temperature
    ]
    compose(ODESystem(eqs, t, vars, para;name),port)
end
function Compressor(;name)
	@named inport = Port()
    @named outport = Port()
	para = @parameters begin
		p1 = 101020, 		[description = "Ambient Pressure (Pa)"]
		T1 = 300, 			[description = "Ambient Temperature (K)"]
		m_cas_comp = 105,   [description = "Mass Casing Compressor (kg)"]
		c_casing = 500, 	[description = "Specific Heat Casing Compressor (J/kg·K)"]
		ṁ_comp = 0.4, 		[description = "Mass Flow Rate Compressor (kg/s)"]
		η_c_is = 0.85, 		[description = "Isentropic Efficiency"]
		π_c = 2.0, 			[description = "Compressor Pressure Ratio"]
		κa = cp(inport.T)/cv(inport.T)
		γ_a = (κa-1) / κa
	end
	vars = @variables begin
		T2(t), 				[description = "Temperature outlet compressor (K)"]
	end
	
	eqs = [ D(T2) ~ (1 / (m_cas_comp * c_casing)) * (ṁ_comp *(h(outport.p,outport.T) 		- h(outport.p, T2)))
			outport.p ~ π_c * inport.p
			outport.T ~ inport.T * (1 + (1 / η_c_is) * (π_c^(γ_a) - 1))
            outport.mdot ~ inport.mdot
			outport.P ~ inport.P - (ṁ_comp * (h(outport.p, outport.T) - h(inport.p, inport.T)))]    
   compose(ODESystem(eqs, t, vars, para;name), inport, outport)
end
function Recuperator(; name, As = 3)
    @named inport_hot = Port()
    @named inport_cold = Port()
    @named outport_hot = Port()
    @named outport_cold = Port()
    
    para = @parameters begin
        A = 10.0 / As,          [description = "Heat transfer area per section (m²)"]
        cp_h = 1000.0,          [description = "Hot fluid heat capacity (J/(kg·K)"]
        cp_c = 1000.0,          [description = "Cold fluid heat capacity (J/(kg·K)"]
        m_dot_h = 0.65,         [description = "Hot side mass flow rate (kg/s)"]
        m_dot_c = 0.65,         [description = "Cold side mass flow rate (kg/s)"]
        h_h = 100.0,            [description = "Hot side heat transfer coefficient (W/m²K)"]
        h_c = 100.0,            [description = "Cold side heat transfer coefficient (W/m²K)"]
        m_w = 185.0,            [description = "Wall mass per section (kg)"]
        c_w = 500.0,            [description = "Wall heat capacity (J/(kg·K)"]
        σ_h = 0.98,             [description = "Hot side pressure loss coefficient"]
        σ_c = 0.95,             [description = "Cold side pressure loss coefficient"]
    end
    
    vars = @variables begin
        T_h_out(t)[1:As],       [description = "Hot fluid outlet temperatures (K)"]
        T_c_out(t)[1:As],       [description = "Cold fluid outlet temperatures (K)"]
        T_w(t)[1:As],           [description = "Wall temperatures (K)"]
    end

    eqs = []
    
    # Temperature equations for each section
    for i in 1:As
		if i == 1
			T_h_in_i = inport_hot.T
		else
			T_h_in_i = T_h_out[i-1]
		end
		
		if i == As
			T_c_in_i = inport_cold.T
		else
			T_c_in_i = T_c_out[i+1]
		end
        # Hot fluid energy balance
        push!(eqs, T_h_out[i] ~ (m_dot_h * cp_h * T_h_in_i - 0.5*h_h*A*T_h_in_i + h_h*A*T_w[i]) / (m_dot_h * cp_h + 0.5*h_h*A))

        # Cold fluid energy balance (counter-flow)
        push!(eqs, T_c_out[i] ~ (m_dot_c * cp_c * T_c_in_i - 0.5*h_c*A*T_c_in_i + h_c*A*T_w[i]) / (m_dot_c * cp_c + 0.5*h_c*A))

        # Wall energy balance (dynamic equation)
        Q_h = m_dot_h * cp_h * (T_h_in_i - T_h_out[i])
        Q_c = m_dot_c * cp_c * (T_c_out[i] - T_c_in_i)
        push!(eqs, D(T_w[i]) ~ (Q_h - Q_c) / (m_w * c_w))
    end

    # System boundary conditions
    push!(eqs,
        outport_hot.T ~ T_h_out[end],   # Connect last hot section to outlet
        outport_cold.T ~ T_c_out[1],    # Connect first cold section to outlet
        outport_hot.p ~ inport_hot.p * σ_h,  # Pressure loss
        outport_cold.p ~ inport_cold.p * σ_c,
        outport_hot.mdot ~ -m_dot_h,    # Mass flow convention
        outport_cold.mdot ~ -m_dot_c,
        inport_hot.mdot ~ m_dot_h,
        inport_cold.mdot ~ m_dot_c)

    compose(ODESystem(eqs, t, vars, para;name), [inport_hot, inport_cold], [outport_hot, outport_cold])
end
function CombustionChamber(;name)
    @named inport = Port()
    @named outport = Port()

    pars = @parameters begin
        η_comb = 0.97,        [description = "Combustor efficiency"]
        σ_comb = 0.945,       [description = "Artificial Pressure Losses"]
        V_comb = 5.5422e-3,   [description = "Combustor Volume (m^3)"]
        LHV = 40.564e6,       [description = "Lower Heating Value (J/kg)"]
        m_cas_comb = 10.0,    [description = "Combustor casing mass (kg)"]
        c_casing = 0,         [description = "Specific heat capacity of casing (J/kg·K)"]
        Tf = 303.15,          [description = "Fuel Temperature (K)"]
        pf = 600000,          [description = "Fuel Pressure (Pa) (6 bar)"]
        ρf = 0.8217,          [description = "Fuel Density (kg/m^3)"]
        ρ3 = 1.15,            [description = "Density in combustor (kg/m^3)"]
        R = 287,              [description = "Specific gas constant for air (J/kg·K)"]
		Rmol = 8.314
		n = 3
		ṁ_comp = 0.4 
		ṁ_turb = 0.46
		ṁ_fuel = 0.06
    end

    vars = @variables begin
        ρ3(t), [description = "Density in combustor"]
        u3(t), [description = "Internal energy in combustor"]
    end

    # Define equations (steady-state: time derivatives are zero)
    eqs = [ D(ρ3) ~ (ṁ_comp + ṁ_fuel - ṁ_turb) / V_comb,
			D(u3) ~ (1 / (ρ3 * V_comb)) * (ṁ_comp * h(inport.p, inport.T) + ṁ_fuel * hf(pf, Tf) + (η_comb * ṁ_fuel * LHV) - ṁ_turb * h(outport.p, outport.T) - u3*(ṁ_comp + ṁ_fuel - ṁ_turb)),
			u3 ~ u(outport.p, outport.T),
        	outport.p ~ ρ3 * R * outport.T,
        	outport.mdot ~ inport.mdot,
        	outport.P ~ 0]
	
    compose(ODESystem(eqs, t, name=name), inport, outport)
end
function Turbine(;name)
	@named inport = Port()
    @named outport = Port()
	para = @parameters begin
		m_cas_turb = 105, 		[description = "Turbine casing mass (kg)"]
        c_casing = 600, 		[description = "Specific Heat Casing Turbine (J/kg·K)"]
        ṁ_t = 0.46,     		[description = "Mass Flow Turbine (kg/s)"]
        η_turb_is = 0.95, 		[description = "Turbine Isentropic Efficiency"]
        π_t = 2.0, 				[description = "Turbine Pressure Ratio"]
		κg = cp(inport.T)/cv(inport.T)
		γ_g = (κg - 1) / κg
	end
	vars = @variables begin
		T4(t), 					[description = "Turbine Outlet Temperature (K)"]
	end
	
	eqs = [ D(T4) ~ (1 / (m_cas_turb * c_casing)) * (ṁ_t * (h(outport.p, outport.T) 		- h(outport.p, T4)))
			outport.p ~ π_t * inport.p
			outport.T ~ inport.T * (1 - η_turb_is * (1 - π_t^(-γ_g)))
            outport.mdot ~ inport.mdot
			outport.P ~ inport.P + (ṁ_t * (h(inport.p, inport.T) - h(outport.p, outport.T)))]	
   compose(ODESystem(eqs, t, vars, para;name), inport, outport)
end
function RotationalSystem(;name)
	@named inport = Port()
    @named outport = Port()
	para = @parameters begin
        I =  8.3e-3, 			[description = "Inertia (kg.m^2)"]
        P_lb = 1500, 			[description = "Power Bearing Losses (W)"]
        P_le = 10000, 			[description = "Power Electrical Losses (W)"]
		Pload = 70, 			[description = "Power Load (W)"]
	end
	vars = @variables begin

	end
	
	eqs = [ outport.p ~ inport.p
			outport.T ~ inport.T
            outport.mdot ~ inport.mdot
			outport.P ~ inport.P - P_lb - P_le]	
   compose(ODESystem(eqs, t, vars, para;name), inport, outport)
end
function Sink(;name) 
    @named    port = Port()
    para = @parameters begin
        
    end
    vars = @variables begin
        p(t)
		T(t)
        mdot(t) 
		P(t)
     end

   eqs = [
        mdot ~ port.mdot
        p ~ port.p
		T ~ port.T
	    P ~ port.P
   ]
   compose(ODESystem(eqs, t, vars, para;name),port)
end

@named source = Source(source_temperature = start_T, source_pressure = start_p,source_mdot = mdot, source_power = start_P)
@named compressor = Compressor()
@named recuperator = Recuperator()
@named combustionchamber = CombustionChamber()
@named turbine = Turbine()
@named sink = Sink()
eqs = [
    connect(source.port,compressor.inport)
	connect(compressor.outport,turbine.inport)
    connect(turbine.outport,sink.port)
]
@named model = ODESystem(eqs, t, systems=[source,compressor,turbine,sink])

However, I’m facing an issue:

In CombustionChamber model, Extracting Output Temperature (output.T):
I need to extract the output temperature T from the combustion chamber using CoolProp, but the problem is that my DAE system requires this T as an input for another equation. This creates a dependency issue, and I’m not sure how to properly structure my equations to resolve it.
In Recuperator model, it causes this error message :

MethodError: no method matching ModelingToolkit.ODESystem(::Vector{Any}, ::Symbolics.Num, ::Vector{Symbolics.Arr{Symbolics.Num, 1}}, ::Vector{Symbolics.Num}; name::Symbol)

The type `ModelingToolkit.ODESystem` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:

  ModelingToolkit.ODESystem(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any; ...)

   @ ModelingToolkit C:\Users\Andhika\.julia\packages\ModelingToolkit\eiNg3\src\systems\diffeqs\odesystem.jl:184

  ModelingToolkit.ODESystem(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any; ...)

   @ ModelingToolkit C:\Users\Andhika\.julia\packages\ModelingToolkit\eiNg3\src\systems\diffeqs\odesystem.jl:184

  ModelingToolkit.ODESystem(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any; ...)"

I would really appreciate any guidance on how to handle these issues. If needed, I can provide more details or share parts of my code. Thank you in advance!

ufechner7 · February 22, 2025, 3:27pm

Do I understand correctly that this problem is coming from the use of MTK components? I always write my DAE equations without using components and this was never an issue, the symbolic engine or the solver resolves any algebraic loops nicely.

AndhikaTama · February 22, 2025, 3:42pm

Yes, when I ran the simulation directly without using components, the code executed successfully without any errors. However, when I modified the code using ModelingToolkit, the issues appeared. Essentially, I aim to build a system similar to HVAC in JuliaSim, as shown in this example: Tube-Fin HX.

ChrisRackauckas · March 2, 2025, 11:12pm

Likely this is causing it to not hit the simplified constructor.

Can you use @mtkmodel instead?

Topic		Replies	Views
Acausal Modelling and Control Sytem Modelling & Simulations question , control , modelling , controlsystems	11	149	March 18, 2025
ModelingToolkit: Index reduction and linear algebraic loops Modelling & Simulations modelingtoolkit	6	993	October 7, 2021
ModelingToolkit: Building a DAEProblem and solving it with IDA with KLU & Jacobian usage Modelling & Simulations modelingtoolkit	3	988	March 20, 2022
Optimization and Modelingtoolkit Modelling & Simulations	7	1088	July 14, 2021
Using ModelingToolkit.jl to model a matrix ODE Modelling & Simulations question , diffeq , modelingtoolkit	3	965	September 17, 2022

ModelingToolkit.jl: Handling Output/Input Dependency in a DAE Model

Related topics