Acausal Modelling and Control Sytem

AndhikaTama · March 16, 2025, 9:50am

Hello everyone,

I am a new user of Julia ModelingToolkit and ControlSystem. I am currently working on a project with the topic Nonlinear Dynamic Simulation of a MicroGas Turbine Integrated with Concentrated Solar Thermal.

As a first step in building this system, I need to validate my MicroGas Turbine (MGT) model by comparing it with the T100 MGT model. I am performing an Acausal Modeling simulation in Julia ModelingToolkit. However, I am facing challenges in implementing the controller for the MGT simulation, which follows the following scheme:

The modular code I have written so far is as follows, and it only lacks the controller block:

using ControlSystems,PlutoUI
using ModelingToolkit, DifferentialEquations, CoolProp, Plots
using ModelingToolkitStandardLibrary
using ModelingToolkitStandardLibrary.Blocks

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

PropsSI2(out::AbstractString, fluid::AbstractString) = CoolProp.PropsSI(out, fluid)
@register_symbolic PropsSI2(out::AbstractString, 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

function R_air() 
	R_air() = PropsSI2("GAS_CONSTANT", fluid)
	return PropsSI2("GAS_CONSTANT", 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)
@register_symbolic R_air()



@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

@connector  function Shaft(;name) 
vars = @variables begin 
    ω(t)
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, ṁ_comp = 0.5, η_c_is = 0.85, π_c = 3.0)
	@named inport = Port()
    @named outport = Port()
	para = @parameters begin
		m_cas_comp = 105,   [description = "Mass Casing Compressor (kg)"]
		c_casing = 500, 	[description = "Specific Heat Casing Compressor (J/kg·K)"]
		κ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, A_total = 10000.0, N = 3)
    @named inport_hot = Port()
    @named inport_cold = Port()
    @named outport_hot = Port()
    @named outport_cold = Port()
    
    # Parameter sistem didefinisikan sebagai NamedTuple (nilai numerik)
    para = (
        cp_h = 500.0,
        cp_c = 500.0,
        m_dot_h = 0.508694,
        m_dot_c = 0.5,
        h_h = 100.0,
        h_c = 100.0,
        m_w = 185.0,
        c_w = 500.0,
        σ_h = 0.98,
        σ_c = 0.95,
        A = A_total / N
    )

    @independent_variables t
    @variables T_w(t)[1:3]
    @variables T_h_out(t)[1:3]
    @variables T_c_out(t)[1:3]
    
    # Inisialisasi vektor persamaan dengan tipe Equation
    eqs = Equation[]
    
    for i in 1:3
        # Definisikan suhu masuk fluida
        T_h_in_i = (i == 1) ? inport_hot.T : T_h_out[i-1]
        T_c_in_i = (i == 3) ? inport_cold.T : T_c_out[i+1]
        
        # Persamaan untuk fluida panas
        push!(eqs, T_h_out[i] ~ ((para.m_dot_h * para.cp_h * T_h_in_i - 0.5 * para.h_h * para.A) * T_h_in_i + para.h_h * para.A * T_w[i]) /
                                ((para.m_dot_h * para.cp_h * T_h_in_i) + (0.5 * para.h_h * para.A)))
        # Persamaan untuk fluida dingin
        push!(eqs, T_c_out[i] ~ ((para.m_dot_c * para.cp_c * T_c_in_i - 0.5 * para.h_c * para.A) * T_c_in_i + para.h_c * para.A * T_w[i]) /
                                ((para.m_dot_c * para.cp_c * T_c_in_i) + (0.5 * para.h_c * para.A)))
        
        # Persamaan diferensial untuk perubahan suhu dinding
		Q_h = para.m_dot_h * (h(inport_hot.p, T_h_in_i) - h(inport_hot.p, T_h_out[i]))
		Q_c = para.m_dot_c * (h(inport_cold.p, T_c_out[i]) - h(inport_cold.p, T_c_in_i))
        push!(eqs, Differential(t)(T_w[i]) ~ (Q_h - Q_c) / ((para.m_w / N) * para.c_w))
    end
    
    # Kondisi batas port
    push!(eqs, outport_hot.T ~ T_h_out[3])
    push!(eqs, outport_cold.T ~ T_c_out[1])
    push!(eqs, outport_hot.p ~ inport_hot.p * para.σ_h)
    push!(eqs, outport_cold.p ~ inport_cold.p * para.σ_c)
    push!(eqs, outport_hot.mdot ~ -para.m_dot_h)
    push!(eqs, outport_cold.mdot ~ -para.m_dot_c)
    push!(eqs, outport_hot.P ~ inport_hot.P)
    push!(eqs, outport_cold.P ~ inport_cold.P)
    
    # Gabungkan semua variabel keadaan ke dalam satu vektor
    flat_vars = vcat(T_w..., T_h_out..., T_c_out...)
    
    # Buat sistem ODE
    sys = ODESystem(eqs, t, flat_vars, []; name=name)
    
    # Komposisi sistem dengan semua port dalam satu vektor
    ports = [inport_hot, inport_cold, outport_hot, outport_cold]
    comp_sys = compose(sys, ports)
    
    return comp_sys
end

function CombustionChamber(;name)
    @named inport = Port()
    @named outport = Port()
	@named ṁ_fuel = RealInput()

    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 = 500,         [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)"]
		ṁ_comp = 0.5 
		ṁ_turb = 0.50627
		xṁ_fuel(t) = 0.00627
    end

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

    # Define equations (steady-state: time derivatives are zero)
    eqs = [ D(ρ3) ~ (ṁ_comp + xṁ_fuel - ṁ_turb) / V_comb,
			D(T3) ~ (1 / (ρ3 * V_comb*cv(outport.T))) * (ṁ_comp * h(inport.p, inport.T) + xṁ_fuel * hf(pf, Tf) + (η_comb * xṁ_fuel * LHV) - ṁ_turb * h(outport.p, outport.T) - (cv(T3)*T3)*(ṁ_comp + xṁ_fuel - ṁ_turb)),
			D(TIT) ~ (ṁ_turb * (h(outport.p, outport.T) - h(outport.p, TIT))) / (m_cas_comb * c_casing),
        	outport.p ~ ρ3 * R_air() * outport.T,
        	outport.mdot ~ inport.mdot,
        	outport.P ~ inport.P,
			outport.T ~ T3,
			xṁ_fuel ~ ṁ_fuel.u]
	
    compose(ODESystem(eqs, t, name=name), inport, outport)
end

function Turbine(;name, ṁ_t = 0.508694, η_turb_is = 0.95, π_t = 3.4)
	@named inport = Port()
    @named outport = Port()
	@named inshaft = Shaft()
	@named outshaft = Shaft()
	@named T4 = RealOutput()
	
	para = @parameters begin
		m_cas_turb = 122.2, 		[description = "Turbine casing mass (kg)"]
        c_casing = 500, 		[description = "Specific Heat Casing Turbine (J/kg·K)"]
		κg = cp(inport.T)/cv(inport.T)
		γ_g = (κg - 1) / κg
	end
	
	vars = @variables begin
		xT4(t), 					[description = "Turbine Outlet Temperature (K)"]
	end
	
	eqs = [ D(xT4) ~ (1 / (m_cas_turb * c_casing)) * (ṁ_t * (h(outport.p, outport.T) 		- h(outport.p, xT4)))
			outport.p ~ inport.p / π_t 
			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)))
			outshaft.ω ~ inshaft.ω
			xT4 ~ T4.u]	

   compose(ODESystem(eqs, t, vars, para;name), inport, outport, inshaft, outshaft)
end

function RotationalSystem(;name)
	@named inport = Port()
    @named outport = Port()
    @named outshaft = Shaft()
	@named Pload = RealInput()
	
	para = @parameters begin
        I =  8.3e-3, 			[description = "Inertia (kg.m^2)"]
		xPload(t) = 70000
        P_lb = 21.33e-3 * (outshaft.ω * (60/(2*pi))),			[description = "Power Bearing Losses (W)"]
        P_le = -5.0802e-3 * (xPload/1000)^3 + 2.3891 * (xPload/1000)^2 − 2.3492e2 * (xPload/1000) + 8.3781e3            , 			[description = "Power Electrical Losses (W)"]
	end
	
	vars = @variables begin
		ω(t)
	end
	
	eqs = [ D(ω) ~ (1 / (I * ω)) * (outport.P/1000)
			outport.p ~ inport.p
			outport.T ~ inport.T
            outport.mdot ~ inport.mdot
			outport.P ~ inport.P - xPload - P_lb - P_le
			outshaft.ω ~ ω
			xPload ~ Pload.u]	
   compose(ODESystem(eqs, t, vars, para;name), inport, outport, outshaft)
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

start_T = 300
start_p = 101325
mdot = 0
start_P = 0

function MicroGasTurbine(;name = :MicroGasTurbine)
	@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 shaft = RotationalSystem()
	@named sink = Sink()
	
	eqs = [ connect(source.port,compressor.inport)
			connect(compressor.outport, recuperator.inport_cold)
			connect(recuperator.outport_cold, combustionchamber.inport)
			connect(combustionchamber.outport, turbine.inport)
			connect(turbine.outport, shaft.inport)
			connect(turbine.outport, recuperator.inport_hot)
			connect(recuperator.outport_hot, sink.port)
			connect(shaft.outshaft, turbine.inshaft)]
		
   	return compose(ODESystem(eqs, t, systems = [source,compressor,combustionchamber,recuperator,turbine, shaft, sink], name = name))
	
end

@named mgt = MicroGasTurbine()

let
	@named mgt = MicroGasTurbine()
	ssys = structural_simplify(mgt)
	u0 = [compressor.T2 => 300, combustionchamber.ρ3 => 10,combustionchamber.T3 => 300, combustionchamber.TIT => 300, recuperator.T_w[1] => 157, recuperator.T_w[2] => 403.5, recuperator.T_w[3] => 650, turbine.xT4 => 300, shaft.ω => 5850]
	tspan = (0.0, 1000.0)
	prob = ODEProblem(ssys, u0, tspan)
	sol = solve(prob, saveat = 1.0)
	plot(sol, layout = 9)
end

Is there anyone here who can provide guidance on designing the controller for this system? I would really appreciate your help!

Thank you in advance!

baggepinnen · March 16, 2025, 9:57am

I’m afraid your ask is too general to get much help. Do you have any concrete questions you’d like answered? Have you tried anything and failed?

AndhikaTama · March 16, 2025, 10:12am

Yes, I watched your video
[link: https://youtu.be/l3odTEy69UU?si=IPuIH_raOwvxh1C-].

I am trying to implement a P controller (as a test) and integrate it into the system as follows:

rṁ_fuel = 0.008
function MicroGasTurbineWithP(;name = :MicroGasTurbineWithP)
	@named ref = Constant(k=650)
	@named controller = Gain(k=20)
	@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 shaft = RotationalSystem()
	@named sink = Sink()
	@named feedback = Feedback()
	
	eqs = [ connect(source.port,compressor.inport)
			connect(compressor.outport, recuperator.inport_cold)
			connect(recuperator.outport_cold, combustionchamber.inport)
			connect(combustionchamber.outport, turbine.inport)
			connect(turbine.outport, shaft.inport)
			connect(turbine.outport, recuperator.inport_hot)
			connect(recuperator.outport_hot, sink.port)
			connect(shaft.outshaft, turbine.inshaft)
			connect(ref.output, feedback.input1)
			connect(turbine.xT4, feedback.input2)
			connect(feedback.output, controller.input)
			connect(controller.output, combustionchamber.xṁ_fuel)]
		
   	return compose(ODESystem(eqs, t, systems = [source,compressor,combustionchamber,recuperator,turbine, shaft, sink, ref, controller, feedback], name = name))
	
end

let
	@named model = MicroGasTurbineWithP()
	ssys = structural_simplify(model)
	tspan = (0.0, 1000.0)
	prob = ODEProblem(ssys, [], tspan)
	sol = solve(prob, Rodas5p())
	plot(sol, layout=2)
	hline!([rṁ_fuel prob[complete(model).ref.k]], label="ref")
end

However, when I run the simulation, I encounter an error, and the following message appears:

Error message from SymbolicUtils
None Sym BasicSymbolic doesn't have a name

Stack trace
Here is what happened, the most recent locations are first:

error(s::String)
from 
julia → error.jl:35
 
nameof(s::SymbolicUtils.BasicSymbolic{Real})
from 
SymbolicUtils → types.jl:271
nameof(n::Symbolics.Num)
from 
Symbolics → num.jl:97
#unique#428(f::typeof(nameof), C::Tuple{…}; seen::Nothing) ...show types...
from 
julia → set.jl:315
 
unique
from 
set.jl:301
connect(::Symbolics.Num, ::ModelingToolkit.ODESystem)
from 
Symbolics → equations.jl:18
MicroGasTurbineWithP(; name::Symbol)
from 
Other cell: line 13
	@named feedback = Feedback()
	
	eqs = [ connect(source.port,compressor.inport)
			connect(compressor.outport, recuperator.inport_cold)
			connect(recuperator.outport_cold, combustionchamber.inport)
Show more...

Is there anything wrong with my code?

baggepinnen · March 16, 2025, 10:35am

Your example cannot be run here, can you please provide a full running example that reproduces the error? I’t very hard to find errors like this by only looking at the code.

AndhikaTama · March 16, 2025, 2:11pm

Sure the full code is this

using ControlSystems,PlutoUI
using ModelingToolkit, DifferentialEquations, CoolProp, Plots
using ModelingToolkitStandardLibrary
using ModelingToolkitStandardLibrary.Blocks

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

PropsSI2(out::AbstractString, fluid::AbstractString) = CoolProp.PropsSI(out, fluid)
@register_symbolic PropsSI2(out::AbstractString, 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

function R_air() 
	R_air() = PropsSI2("GAS_CONSTANT", fluid)
	return PropsSI2("GAS_CONSTANT", 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)
@register_symbolic R_air()



@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

@connector  function Shaft(;name) 
vars = @variables begin 
    ω(t)
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, ṁ_comp = 0.5, η_c_is = 0.85, π_c = 3.0)
	@named inport = Port()
    @named outport = Port()
	para = @parameters begin
		m_cas_comp = 105,   [description = "Mass Casing Compressor (kg)"]
		c_casing = 500, 	[description = "Specific Heat Casing Compressor (J/kg·K)"]
		κ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, A_total = 10000.0, N = 3)
    @named inport_hot = Port()
    @named inport_cold = Port()
    @named outport_hot = Port()
    @named outport_cold = Port()
    
    # Parameter sistem didefinisikan sebagai NamedTuple (nilai numerik)
    para = (
        cp_h = 500.0,
        cp_c = 500.0,
        m_dot_h = 0.508694,
        m_dot_c = 0.5,
        h_h = 100.0,
        h_c = 100.0,
        m_w = 185.0,
        c_w = 500.0,
        σ_h = 0.98,
        σ_c = 0.95,
        A = A_total / N
    )

    @independent_variables t
    @variables T_w(t)[1:3]
    @variables T_h_out(t)[1:3]
    @variables T_c_out(t)[1:3]
    
    # Inisialisasi vektor persamaan dengan tipe Equation
    eqs = Equation[]
    
    for i in 1:3
        # Definisikan suhu masuk fluida
        T_h_in_i = (i == 1) ? inport_hot.T : T_h_out[i-1]
        T_c_in_i = (i == 3) ? inport_cold.T : T_c_out[i+1]
        
        # Persamaan untuk fluida panas
        push!(eqs, T_h_out[i] ~ ((para.m_dot_h * para.cp_h * T_h_in_i - 0.5 * para.h_h * para.A) * T_h_in_i + para.h_h * para.A * T_w[i]) /
                                ((para.m_dot_h * para.cp_h * T_h_in_i) + (0.5 * para.h_h * para.A)))
        # Persamaan untuk fluida dingin
        push!(eqs, T_c_out[i] ~ ((para.m_dot_c * para.cp_c * T_c_in_i - 0.5 * para.h_c * para.A) * T_c_in_i + para.h_c * para.A * T_w[i]) /
                                ((para.m_dot_c * para.cp_c * T_c_in_i) + (0.5 * para.h_c * para.A)))
        
        # Persamaan diferensial untuk perubahan suhu dinding
		Q_h = para.m_dot_h * (h(inport_hot.p, T_h_in_i) - h(inport_hot.p, T_h_out[i]))
		Q_c = para.m_dot_c * (h(inport_cold.p, T_c_out[i]) - h(inport_cold.p, T_c_in_i))
        push!(eqs, Differential(t)(T_w[i]) ~ (Q_h - Q_c) / ((para.m_w / N) * para.c_w))
    end
    
    # Kondisi batas port
    push!(eqs, outport_hot.T ~ T_h_out[3])
    push!(eqs, outport_cold.T ~ T_c_out[1])
    push!(eqs, outport_hot.p ~ inport_hot.p * para.σ_h)
    push!(eqs, outport_cold.p ~ inport_cold.p * para.σ_c)
    push!(eqs, outport_hot.mdot ~ -para.m_dot_h)
    push!(eqs, outport_cold.mdot ~ -para.m_dot_c)
    push!(eqs, outport_hot.P ~ inport_hot.P)
    push!(eqs, outport_cold.P ~ inport_cold.P)
    
    # Gabungkan semua variabel keadaan ke dalam satu vektor
    flat_vars = vcat(T_w..., T_h_out..., T_c_out...)
    
    # Buat sistem ODE
    sys = ODESystem(eqs, t, flat_vars, []; name=name)
    
    # Komposisi sistem dengan semua port dalam satu vektor
    ports = [inport_hot, inport_cold, outport_hot, outport_cold]
    comp_sys = compose(sys, ports)
    
    return comp_sys
end

function CombustionChamber(;name)
    @named inport = Port()
    @named outport = Port()
	@named ṁ_fuel = RealInput()

    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 = 500,         [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)"]
		ṁ_comp = 0.5 
		ṁ_turb = 0.50627
		xṁ_fuel(t) = 0.00627
    end

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

    # Define equations (steady-state: time derivatives are zero)
    eqs = [ D(ρ3) ~ (ṁ_comp + xṁ_fuel - ṁ_turb) / V_comb,
			D(T3) ~ (1 / (ρ3 * V_comb*cv(outport.T))) * (ṁ_comp * h(inport.p, inport.T) + xṁ_fuel * hf(pf, Tf) + (η_comb * xṁ_fuel * LHV) - ṁ_turb * h(outport.p, outport.T) - (cv(T3)*T3)*(ṁ_comp + xṁ_fuel - ṁ_turb)),
			D(TIT) ~ (ṁ_turb * (h(outport.p, outport.T) - h(outport.p, TIT))) / (m_cas_comb * c_casing),
        	outport.p ~ ρ3 * R_air() * outport.T,
        	outport.mdot ~ inport.mdot,
        	outport.P ~ inport.P,
			outport.T ~ T3,
			xṁ_fuel ~ ṁ_fuel.u]
	
    compose(ODESystem(eqs, t, name=name), inport, outport)
end

function Turbine(;name, ṁ_t = 0.508694, η_turb_is = 0.95, π_t = 3.4)
	@named inport = Port()
    @named outport = Port()
	@named inshaft = Shaft()
	@named outshaft = Shaft()
	@named T4 = RealOutput()
	
	para = @parameters begin
		m_cas_turb = 122.2, 		[description = "Turbine casing mass (kg)"]
        c_casing = 500, 		[description = "Specific Heat Casing Turbine (J/kg·K)"]
		κg = cp(inport.T)/cv(inport.T)
		γ_g = (κg - 1) / κg
	end
	
	vars = @variables begin
		xT4(t), 					[description = "Turbine Outlet Temperature (K)"]
	end
	
	eqs = [ D(xT4) ~ (1 / (m_cas_turb * c_casing)) * (ṁ_t * (h(outport.p, outport.T) 		- h(outport.p, xT4)))
			outport.p ~ inport.p / π_t 
			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)))
			outshaft.ω ~ inshaft.ω
			xT4 ~ T4.u]	

   compose(ODESystem(eqs, t, vars, para;name), inport, outport, inshaft, outshaft)
end

function RotationalSystem(;name)
	@named inport = Port()
    @named outport = Port()
    @named outshaft = Shaft()
	@named Pload = RealInput()
	
	para = @parameters begin
        I =  8.3e-3, 			[description = "Inertia (kg.m^2)"]
		xPload(t) = 70000
        P_lb = 21.33e-3 * (outshaft.ω * (60/(2*pi))),			[description = "Power Bearing Losses (W)"]
        P_le = -5.0802e-3 * (xPload/1000)^3 + 2.3891 * (xPload/1000)^2 − 2.3492e2 * (xPload/1000) + 8.3781e3            , 			[description = "Power Electrical Losses (W)"]
	end
	
	vars = @variables begin
		ω(t)
	end
	
	eqs = [ D(ω) ~ (1 / (I * ω)) * (outport.P/1000)
			outport.p ~ inport.p
			outport.T ~ inport.T
            outport.mdot ~ inport.mdot
			outport.P ~ inport.P - xPload - P_lb - P_le
			outshaft.ω ~ ω
			xPload ~ Pload.u]	
   compose(ODESystem(eqs, t, vars, para;name), inport, outport, outshaft)
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

start_T = 300
start_p = 101325
mdot = 0
start_P = 0

function MicroGasTurbineWithP(;name = :MicroGasTurbineWithP)
	@named ref = Constant(k=650)
	@named controller = Gain(k=20)
	@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 shaft = RotationalSystem()
	@named sink = Sink()
	@named feedback = Feedback()
	
	eqs = [ connect(source.port,compressor.inport)
			connect(compressor.outport, recuperator.inport_cold)
			connect(recuperator.outport_cold, combustionchamber.inport)
			connect(combustionchamber.outport, turbine.inport)
			connect(turbine.outport, shaft.inport)
			connect(turbine.outport, recuperator.inport_hot)
			connect(recuperator.outport_hot, sink.port)
			connect(shaft.outshaft, turbine.inshaft)
			connect(ref.output, feedback.input1)
			connect(turbine.xT4, feedback.input2)
			connect(feedback.output, controller.input)
			connect(controller.output, combustionchamber.xṁ_fuel)]
		
   	return compose(ODESystem(eqs, t, systems = [source,compressor,combustionchamber,recuperator,turbine, shaft, sink, ref, controller, feedback], name = name))
	
end

rṁ_fuel = 0.008

let
	@named model = MicroGasTurbineWithP()
	ssys = structural_simplify(model)
	tspan = (0.0, 1000.0)
	prob = ODEProblem(ssys, [], tspan)
	sol = solve(prob, Rodas5p())
	plot(sol, layout=2)
	hline!([rṁ_fuel prob[complete(model).ref.k]], label="ref")
end

And it produces an error like in this message warning

baggepinnen · March 16, 2025, 4:17pm

It looks like you’re trying to connect a variable to a connector, the variable is the first argument of the call to connect. This is likely the offending variable e
xT4

AndhikaTama · March 16, 2025, 4:27pm

How can I resolve this issue if I want to regulate the turbine model’s T4 output using a reference temperature by controlling the fuel mass flow rate (ṁ_fuel) from the combustion chamber model?

baggepinnen · March 16, 2025, 4:29pm

Either create a connector in the turbine component or use an equation instead of a connect statement

@cryptic.ax should connect perhaps handle when one side is a connector and the other is a variable?

cryptic.ax · March 17, 2025, 5:31am

We could. The easiest way is to use the causal variable connect which simply creates an equality and validates that the connector(s) have a single variable.

AndhikaTama · March 18, 2025, 2:51am

After studying it, it seems that the error occurred because I did not use JuliaSimControl. Therefore, can I still design a controller for acausal modeling without using JuliaSimControl? And have you ever made a tutorial on that?

baggepinnen · March 18, 2025, 5:13am

I very much doubt that, your problem was an illegal connection, and JSControl does nothing to make that legal.

Yes, you can design such a controller with pen and paper or any other method of your liking.

Several

The one you linked is relying very little on JSControl
Batch linearization and gain scheduling · ControlSystemsMTK Documentation
Modeling for control using ModelingToolkit Â· JuliaSimControl
DC Motor with Speed Controller · ModelingToolkitStandardLibrary.jl
https://www.youtube.com/playlist?list=PLC0QOsNQS8hZtOQPHdtul3kpQwMOBL8Qc

AndhikaTama · March 18, 2025, 6:12am

Thank you! i’ll check the references and give you an update on the results.

AndhikaTama · April 28, 2025, 10:27am

Hello Fredrik,
I’m writing to update you on my project progress. I’ve successfully resolved the previous issue you mentioned about the illegal connection between variables and have established the correct connections.

However, I’ve encountered a new challenge regarding the Pload variable. As shown in my system diagram :

As you can see, the system need Pload to serve both as:

An output (to provide data from Pload-1 for calculations)
An input (to feed updated Pload values back into the system)

This creates a circular dependency where we need to:
Use the previous timestep’s Pload value as output data
Simultaneously input a new Pload value for the current timestep

Key Question:
Is it possible to implement a variable that can:

Maintain its value from the previous calculation step (as output)
Accept updated input values for the current step
All within the same variable or through a proper architectural pattern?

Here is the key portion of my code (focusing only on the essential parts involving Pload and system integration):

function RotationalSystem(; name)
    @named InPowPortFromComp = PowerPort()
    @named InPowPortFromTurb = PowerPort()
    @named OutPowPort = PowerPort()
    @named outshaft = Shaft()
	@named ω = RealOutput()
	@named Pload = RealInput()

    para = @parameters begin
        I =  8.3e-3                 # Inertia
        P_lb = 21.33e-3 * (outshaft.ω xω)     # Power bearing losses
        P_le = 3473 + 0.0812*(InPowPortFromTurb.P-InPowPortFromComp.P) +
              2.67e-7 * (InPowPortFromTurb.P-InPowPortFromComp.P)^2 -
              4.48e-12 * (InPowPortFromTurb.P-InPowPortFromComp.P)^3   # Electrical losses
    end

    vars = @variables begin
        xω(t) = 7000
        xPload(t) = 100000        
    end

    eqs = [
        D(xω) ~ (1 / (I * xω)) * (OutPowPort.P / 1000),
        outshaft.ω ~ xω,
		OutPowPort.P ~ (InPowPortFromTurb.P - InPowPortFromComp.P) - xPload - P_lb - P_le,

		#System's Sensor
		xω ~ xω.u
		xPload ~ Pload.u,
    ]

    compose(ODESystem(eqs, t, vars, para; name),
        InPowPortFromComp, InPowPortFromTurb, OutPowPort,
        outshaft, ω, Pload)
end

function MicroGasTurbineWithController(; name = :MicroGasTurbineWithController)
    # === Komponen Plant (Sama seperti sebelumnya) ===
    @named sourceT             = SourceTemperature(source_temperature = 300)
    @named sourcep             = SourcePressure(source_pressure = 101320)
    @named compressor          = CompressorWithMap()
    @named turbine             = TurbineWithMap()
    @named recuperator         = Recuperator()
    @named combustionchamber   = CombustionChamber()
    @named shaft               = RotationalSystem()
    @named wasteheat           = WasteHeat()
    @named ambientoutpressure  = AmbientOutPressure()

	# === Komponen Controller ===
	@named pid_power = PID(k=400, Ti=0.5, Td=1, Nd=10, int__x=0, der__x=0)  # PID untuk power control (menghasilkan ΔN)
	@named pid_speed = PID(k=300, Ti=0.3, Td=0.5, Nd=10, int__x=0, der__x=0)  # PID untuk speed control (output ke Pload)
	@named pid_TOT = PID(k=200, Ti=1.0, Td=0.1, Nd=10, int__x=0, der__x=0)  # PID khusus TOT compensation

    # === Sinyal dan Sensor ===
    @named r_power    = Step(start_time = 0, height=10000)  # Permintaan daya (contoh: step 1000 W)
	
    # === Blok Kustom ===
    @named power_limiter = Limiter(;   name=:power_limiter, y_max=100000, y_min=70000)
    @named TOT_ref_logic     = TOT_Reference_Logic()    # Hitung TOT_ref berdasarkan Tamb
    @named integrator_TOT    = Integrator(;name, k = 1, x = 0.0) # Integrasi kompensasi TOT

    @named LUT_Pdem_to_RefSpeed = functionBlock_LUT_Pdem()
	@named LUT_Speed_to_RefFuel = functionBlock_LUT_Fuel()

    # === Blok Adder untuk Power Control ===
	@named sum_power = Add(; name=:sum_power, k1 = 1.0, k2 = 1.0)
    @named sum_fuel1 = Add(; name=:sum_fuel1, k1 = 1.0, k2 = 1.0)
    @named sum_fuel2 = Add(; name=:sum_fuel2, k1 = 1.0, k2 = 1.0)

	# Feedback
	@named feedbackpower = Feedback()
	@named feedbackspeed = Feedback()
	@named feedbackTOT = Feedback()
	
    # ======================
    # Persamaan Koneksi
    # ======================
    eqs = [
        # --- Plant Connections (Tetap sama) ---
        connect(sourceT.TempPort, compressor.InTempPort)
        connect(sourcep.PresPort, compressor.InPresPort)
        connect(shaft.outshaft, compressor.shaftin)

        # Ke Recuperator
        connect(compressor.OutPresPort, recuperator.InPresPort_cold)
        connect(compressor.OutTempPort, recuperator.InTempPort_cold)
        connect(compressor.OutMdotPort, recuperator.InMdotPort_cold)
        connect(turbine.OutPresPort, recuperator.InPresPort_hot)
        connect(turbine.OutMdotPort, recuperator.InMdotPort_hot)
        connect(turbine.OutTempPort, recuperator.InTempPort_hot)

        # Ke Combustion Chamber
        connect(recuperator.OutPresPort_cold, combustionchamber.InPresPort)
        connect(recuperator.OutTempPort_cold, combustionchamber.InTempPort)
        connect(compressor.OutMdotPort, combustionchamber.InMdotPortFromCompressor)
        connect(turbine.OutMdotPort, combustionchamber.InMdotPortFromTurbine)

        # Ke Turbine
        connect(combustionchamber.OutTempPort, turbine.InTempPort)
        connect(combustionchamber.OutPresPort, turbine.InPresPort)
        connect(shaft.outshaft, turbine.shaftin) 
            
        # Ke Shaft
        connect(compressor.OutPowPort, shaft.InPowPortFromComp)
        connect(turbine.OutPowPort, shaft.InPowPortFromTurb)
            
        # Heat Waste
        connect(recuperator.OutTempPort_hot, wasteheat.TempPort)
        connect(recuperator.OutPresPort_hot, ambientoutpressure.PresPort)
		
        # --- Power Control Loop (Two-Stage) ---
        # Langkah 1: Rate limiter pada Pdem
        connect(r_power.output, power_limiter.input)
        # Langkah 2: Lookup table mengubah Pdem menjadi kecepatan referensi awal (N_ref0)
        connect(power_limiter.output, LUT_Pdem_to_RefSpeed.input)
        # Langkah 3: PID Power menghasilkan koreksi ΔN, menggunakan Pdem (setpoint) dan Pload (feedback)
        connect(power_limiter.output, feedbackpower.input1)
        connect(shaft.Pload, feedbackpower.input2)
		connect(feedbackpower.output, pid_power.err_input)
        # Langkah 4: Penjumlahan: LUT output (N_ref0) + koreksi ΔN
        connect(LUT_Pdem_to_RefSpeed.output, sum_power.input1)
        connect(pid_power.ctr_output, sum_power.input2)
        # Langkah 5: Output penjumlahan N_ref dikirim sebagai referensi ke PID Speed
        connect(sum_power.output, feedbackspeed.input1)
        # Langkah 6: PID Speed membandingkan referensi N_ref dengan kecepatan aktual dan mengatur Pload
        connect(shaft.ω, feedbackspeed.input2)
		connect(feedbackspeed.output, pid_speed.err_input)
        connect(pid_speed.ctr_output, shaft.Pload)

        # --- Fuel Control Loop (TOT Compensation) ---
        # Pengolahan sinyal TOT
        connect(sourceT.T_amb_sensor, TOT_ref_logic.Tamb)
        connect(TOT_ref_logic.TOT_ref, feedbackTOT.input1)
        connect(turbine.T4, feedbackTOT.input2)
		connect(feedbackTOT.output, pid_TOT.err_input)
        connect(pid_TOT.ctr_output, integrator_TOT.input)
        # Adder Fuel Stage 1: Gabungkan base fuel dan kompensasi langsung dari PID TOT
        connect(shaft.ω, LUT_Speed_to_RefFuel.input)     # Base fuel dari kecepatan
        connect(LUT_Speed_to_RefFuel.output, sum_fuel1.input1)
        connect(pid_TOT.ctr_output, sum_fuel1.input2)
        # Adder Fuel Stage 2: Gabungkan hasil adder pertama dengan output integrator (aksi integral TOT)
        connect(sum_fuel1.output, sum_fuel2.input1)
        connect(integrator_TOT.output, sum_fuel2.input2)
        # Output akhir fuel control sebagai laju bahan bakar ke combustion chamber
        connect(sum_fuel2.output, combustionchamber.ṁ_fuel)
    ]

    return compose(ODESystem(eqs, t, name=name, systems=[
        sourceT, sourcep, compressor, turbine, recuperator, combustionchamber, shaft,
        wasteheat, ambientoutpressure, pid_power, pid_speed, pid_TOT, r_power, power_limiter, TOT_ref_logic, integrator_TOT, LUT_Pdem_to_RefSpeed, LUT_Speed_to_RefFuel, sum_power, sum_fuel1, sum_fuel2, feedbackpower, feedbackspeed, feedbackTOT
    ]))
end

I’d appreciate your guidance on how to best implement this bidirectional data flow while avoiding algebraic loops or causality issues.
Looking forward to your expert advice!

baggepinnen · April 28, 2025, 10:58am

Pload in this case is what we refer to as a discrete-time variable. The support for discrete-time models in MTK is still limited (very active work in progress), but you cal already now make use of a DiscreteCallback

see particularily periodic events. You would use this by

Adding a time-varying parameter @parameters Pload(t)
Add a periodic event discrete_events = [sample_interval => [Pload ~ f(Pload)]] which is interpreted as Pload(t) = f(Pload(t-1))

It’s not clear to me what parts of your model are supposed to be discrete time and what parts are continuous time, you use continuous-time components, but have a discrete-time input?

AndhikaTama · April 28, 2025, 11:15am

I’m trying to create a system for a microgas turbine following a scheme like this:

Based on this logic, you can see that the process will be continuous-time across all components. However, the logic for power control requires the power from the previous time step to update the Pload applied to the system.
Is using a time-varying parameter approach appropriate for this case?

Matthew1 · April 29, 2025, 9:58am

PID’s should be loop-breakers themselves, no?

Otherwise, one can add a delay: D(P_) ~ (P_ - P)/tau; tau: 1e-6 or something

In Modelica, there is prev(P), but not in MTK I believe.

baggepinnen · April 29, 2025, 10:07am

ModelingToolkit does not yet have the user-facing concept of a discrete-time state variable, and conflates this with a time-varying parameter.

what does the dashed line here signify? Is the entire CS block operating in discrete time? Does this have state variables of its own?

AndhikaTama · May 1, 2025, 5:45am

So basically the CS contains of 2 control system as it shows in the images below.

Here is the explanation from the writer for both controller logic :
3.7.1 Power control
This control loop aims at regulating the mGT speed through a control of the load (Pload)
applied by the generator to the shaft. This load corresponds to the electrical power that is
produced by the mGT.
First, a look-up table links the demanded power (Pdem) to a predefined reference speed.
Afterwards, this reference speed is corrected with a PID regulation of the fault signal between
Pdem and Pload. Finally, the error between the corrected reference speed and the actual shaft
speed is sent to another PID controller that outputs an updated Pload. The power control
scheme is represented in Fig. 3.19.

3.7.2 Fuel control
This control system regulates the amount of fuel that is injected inside the combustion
chamber. It is subdivided into three different loops. The injected fuel rate corresponds to the
least among the three different computed values.
• Acceleration loop: used during the start-up phase to limit the acceleration of the shaft.
• Speed loop: used to ensure that the machine is not over-speeding.
• Temperature loop: used to maintain TOT constant.
However, only the TOT control is implemented in this work, because the two other loops
schemes were not available. According to [28], the TOT loop is used most of the time in normal
conditions.
At first, a compensation fuel requirement is evaluated based on the ambient temperature
(Tamb) and the TOT. A reference value is attributed to TOT depending on Tamb (e.g. 645oC at
15oC). This reference value is selected to maximise the mGT efficiency (i.e. maximization of
temperature maximizes the Carnot efficiency) while respecting the constraint on the maximum
recuperator hot leg inlet temperature. With reduced ambient temperatures, the compression
is achieved more efficiently and less power is consequently required to drive the compressor.
Therefore, the TIT, and thus also the TOT, must be limited to reduce the power generated by
the turbine.
In a second time, a predefined reference fuel requirement is attributed based on the shaft
speed via a look-up table. The TOT compensation fuel is then added to this reference fuel
requirement. This newly defined corrected fuel requirement is finally summed with the integral
of the TOT compensation fuel.
The fuel control scheme is represented in Fig. 3.20. In this work, as mentioned in Section
3.3, the compressor used to increase the natural gas pressure is not modelled. It is therefore
assumed that the input natural gas naturally matches the fuel system requirements prescribed
by Turbec (i.e. 6.0bar [27]).

The writer previously work on this project using Python language with this algorithm logic :

Can I create the same algorithm using the Acausal ModelingToolkit, or should I use other tools in Julia to solve this problem?

ufechner7 · May 1, 2025, 6:03am

With the risk of repeating myself: You can always model your physical system using MTK in continues time and implement your controller using discrete time controllers like DiscretePIDs and couple them by running the phyical system step by step, using the outputs of the integrator as input of your controller and using the output of your controller as parameters of you continues time system.

A (not so simple) example of the step! function can be found here. What you can see in this example is that we run parts of the model less frequently than the main model.

The related MTK functions are defined here.

An example for a controller with an inner and outer loop and a block the linearizes the system can be found here.

I still need to write a tutorial for this topic.

AndhikaTama · May 1, 2025, 6:39am

I can’t really understand what the example is trying to represent. Could you explain in simple terms how the system in the example works?

Topic		Replies	Views
ModelingToolkit.jl: Handling Output/Input Dependency in a DAE Model Modelling & Simulations	3	82	March 2, 2025
ModelingToolkit: Some questions on usage & best practices Modelling & Simulations modelingtoolkit	5	1762	March 20, 2022
Symbols and variables Modelling & Simulations question , package , modelingtoolkit	6	284	October 23, 2023
Acausal modeling and Unreal Engine Modelling & Simulations modelingtoolkit , modelling , games	16	1005	May 8, 2024
MTK linearization: Change of use, or bug? Modelling & Simulations	8	66	August 8, 2024

Acausal Modelling and Control Sytem

Related topics