Acausal Modelling and Control Sytem

[AndhikaTama]

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, ṁ_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)) * (ṁ_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 - (ṁ_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 ṁ_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)"]
		ṁ_comp = 0.5 
		ṁ_turb = 0.50627
		xṁ_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) ~ (ṁ_comp + xṁ_fuel - ṁ_turb) / V_comb,
			D(T3) ~ (1 / (ρ3 * V_comb*cv(outport.T))) * (ṁ_comp * h(inport.p, inport.T) + xṁ_fuel * hf(pf, Tf) + (η_comb * xṁ_fuel * LHV) - ṁ_turb * h(outport.p, outport.T) - (cv(T3)*T3)*(ṁ_comp + xṁ_fuel - ṁ_turb)),
			D(TIT) ~ (ṁ_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,
			xṁ_fuel ~ ṁ_fuel.u]
	
    compose(ODESystem(eqs, t, name=name), inport, outport)
end

function Turbine(;name, ṁ_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)) * (ṁ_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 + (ṁ_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]

@AndhikaTama

Is there anyone here who can provide guidance on designing the controller for this system?

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]

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:

rṁ_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.xṁ_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!([rṁ_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]

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]

@AndhikaTama

(sol, layout = 9)
end

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, ṁ_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)) * (ṁ_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 - (ṁ_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 ṁ_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)"]
		ṁ_comp = 0.5 
		ṁ_turb = 0.50627
		xṁ_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) ~ (ṁ_comp + xṁ_fuel - ṁ_turb) / V_comb,
			D(T3) ~ (1 / (ρ3 * V_comb*cv(outport.T))) * (ṁ_comp * h(inport.p, inport.T) + xṁ_fuel * hf(pf, Tf) + (η_comb * xṁ_fuel * LHV) - ṁ_turb * h(outport.p, outport.T) - (cv(T3)*T3)*(ṁ_comp + xṁ_fuel - ṁ_turb)),
			D(TIT) ~ (ṁ_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,
			xṁ_fuel ~ ṁ_fuel.u]
	
    compose(ODESystem(eqs, t, name=name), inport, outport)
end

function Turbine(;name, ṁ_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)) * (ṁ_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 + (ṁ_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.xṁ_fuel)]
		
   	return compose(ODESystem(eqs, t, systems = [source,compressor,combustionchamber,recuperator,turbine, shaft, sink, ref, controller, feedback], name = name))
	
end

rṁ_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!([rṁ_fuel prob[complete(model).ref.k]], label="ref")
end

And it produces an error like in this message warning

[baggepinnen]

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]

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]

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]

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]

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]

@AndhikaTama

After studying it, it seems that the error occurred because I did not use JuliaSimControl.

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

@AndhikaTama

can I still design a controller for acausal modeling without using JuliaSimControl?

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

@AndhikaTama

have you ever made a tutorial on that?

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]

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