UDE to solve buffer #505

gregorlein · 2024-04-18T12:05:14Z

gregorlein
Apr 18, 2024

Hi,
I'm trying to model a mixing temperature in a buffer with an ude inside my tespy model. The equation is:
t24 = t23 +(t0-t23) * exp(-v23/V)
with v23 being the flow through the buffer, V the volume of the buffer and t0 the initial temperature of the buffer before mixture. t24 is the temperature coming out of the buffer after full mixing. Before solving the model, I don't know the flow v23. So initially I was using a guess for t24 and then iterating to solve the tespy model as well as the equation above. That needs a lot of computational time, so I thought it would make sense to implement it directly within the tespy model? I tried to follow the example in the documentation but I'm not fully understanding it. I get the error 'int' object is not callable in my solve_buffer function... Thanks for any help!

def option_1_mode5_test(q_owd, q_tud, t_supply, t_return_tud, t_return_owd, q_condenser, t_buffer, iterinfo=False): 
    demand_owd = q_owd*-1e6
    demand_tud = q_tud*-1e6
    q_condenser = q_condenser*-1e6

    dap_si = Sink('Sink')
    dap_so = Source('Source')
    me_bypass = Merge('Merge Bypass')
    me_owd = Merge('Merge OWD')
    va_bypass = Valve('Valve Bypass')
    sp_tud = Splitter('Splitter TUD')
    sp_tud_1 = Splitter('Splitter TUD 1')
    va_tud = Valve('Valve TUD')
    cc = CycleCloser('CC')
    pu_tud = Pump('Pump before TUD')
    pu_ev = Pump('Pump before Evaporator')
    he_dap = HeatExchanger('DAP')

    tud = SimpleHeatExchanger('TUD')
    owd = SimpleHeatExchanger('OWD')

    buffer = SimpleHeatExchanger('Buffer')
    sp_buffer = Splitter('Splitter Buffer')
    me_buffer = Merge('Merge Buffer')
    va_buffer = Valve('Valve Buffer')
    
    # Heat Pump
    cd = Condenser('Condenser')
    va_hp = Valve('Expansion Valve')
    cp = Compressor('Compressor')
    ev = HeatExchanger('Evaporator')
    cc_hp = CycleCloser('CC Heat Pump')

    c0 = Connection(dap_so, 'out1', he_dap, 'in1', label='0')
    c1 = Connection(he_dap, 'out1', dap_si, 'in1', label='1')
    c5 = Connection(he_dap, 'out2', me_bypass, 'in1', label='5')
    c7 = Connection(me_bypass, 'out1', cd, 'in2', label='7')
    c8 = Connection(cd, 'out2', cc, 'in1', label='8')
    c9 = Connection(cc, 'out1', pu_tud, 'in1', label='9')
    c10= Connection(pu_tud, 'out1', sp_buffer, 'in1', label='10')
    c12 = Connection(sp_buffer, 'out1', sp_tud, 'in1', label='12')
    c13 = Connection(sp_tud, 'out1', owd, 'in1', label='13')
    c14 = Connection(owd, 'out1', me_owd, 'in1', label='14')
    c15 = Connection(va_bypass, 'out1', me_owd, 'in2', label='15')
    c16 = Connection(sp_tud_1, 'out1', va_bypass, 'in1', label='16')
    c17 = Connection(sp_tud_1, 'out2', me_buffer, 'in1', label='17')
    c18 = Connection(me_buffer, 'out1', me_bypass, 'in2', label='18')    
    c19 = Connection(me_owd, 'out1', pu_ev, 'in1', label='19')
    c20 = Connection(pu_ev, 'out1', ev, 'in1', label='20')
    c21 = Connection(ev, 'out1', he_dap, 'in2', label='21')
    c22 = Connection(sp_buffer, 'out2', va_buffer, 'in1', label='22')
    c23 = Connection(va_buffer, 'out1', buffer, 'in1', label='23')
    c24 = Connection(buffer, 'out1', me_buffer, 'in2', label='24')
    c25 = Connection(sp_tud, 'out2', va_tud, 'in1', label='25')
    c26 = Connection(va_tud, 'out1', tud, 'in1', label='26') 
    c27 = Connection(tud, 'out1', sp_tud_1, 'in1', label='27')

    # HP
    c30 = Connection(cc_hp, 'out1', cp, 'in1', label='30')
    c31 = Connection(cp, 'out1', cd, 'in1', label='31')
    c32 = Connection(cd, 'out1', va_hp, 'in1', label='32')
    c33 = Connection(va_hp, 'out1', ev, 'in2', label='33')
    c34 = Connection(ev, 'out2', cc_hp, 'in1', label='34')

    nw = Network(T_unit='C', p_unit='bar', h_unit='kJ / kg', m_unit = 'kg / s')
    nw.add_conns(c0, c1, c5, c7, c8, c9, c10, c12, c13, c14, c15, c16, c17, c18, c19, c20, c21, c22, c23, c24, c25, c26, c27, c30, c31, c32, c33, c34)

    he_dap.set_attr(pr1=0.98, pr2=0.98)
    tud.set_attr(pr=0.98, Q = demand_tud)
    owd.set_attr(pr=0.98, Q = demand_owd)
    buffer.set_attr(pr=0.98)
    pu_ev.set_attr(eta_s=0.75)
    pu_tud.set_attr(eta_s=0.75)

    # Heat Pump
    cd.set_attr(pr1=0.98, pr2=0.98, Q = q_condenser)
    ev.set_attr(pr1=0.98, pr2=0.98)
    cp.set_attr(eta_s=0.85)

    c0.set_attr(fluid={'water': 1}, T=76.5, p=1, v=375/3600)
    c5.set_attr(fluid={'water': 1}, T=74.5)
    c7.set_attr(p=1)
    c19.set_attr(p=1)
    c8.set_attr(T=t_supply+1)
    c20.set_attr(v=375/3600)
    c27.set_attr(T = t_return_tud+1)
    c14.set_attr(T = t_return_owd+1)

    # Buffer
    #c24.set_attr(T=t_buffer)
    
    # Heat Pump
    # condensation point
    p_cond = CP.PropsSI("P", "Q", 0, "T", t_supply + 5 + 273.15, 'NH3') * 1e-5
    c32.set_attr(p=p_cond)
    # evaporation point
    p_eva = CP.PropsSI("P", "Q", 1, "T", 16 + 273.15, 'NH3') * 1e-5
    c34.set_attr(p=p_eva, x=1, fluid={'NH3':1})
    
    nw.set_attr(iterinfo=iterinfo)

    def solve_buffer(ude): 
        t_tank_initial = ude.params['t_tank_initial']
        volume_tank = ude.params['volume_tank']
        equation = ude.conns[0].calc_T() + (t_tank_initial - ude.conns[0].calc_T()) * np.exp((ude.conns[0].v.val_SI()*3600)/volume_tank) - ude.conns[1].calc_T()
        return equation

    def solve_buffer_deriv(ude):
        c23 = ude.conns[0]
        c24 = ude.conns[1]
        if c23.m.is_var:
            ude.jacobian[c23.m.J_col] = ude.numeric_deriv('m', c23)
        if c23.p.is_var:
            ude.jacobian[c23.p.J_col] = ude.numeric_deriv('p', c23)
        if c23.h.is_var:
            ude.jacobian[c23.h.J_col] = ude.numeric_deriv('h', c23)
        if c24.p.is_var:
            ude.jacobian[c24.p.J_col] = ude.numeric_deriv('p', c24)
        if c24.h.is_var:
            ude.jacobian[c24.h.J_col] = ude.numeric_deriv('h', c24)

    ude = UserDefinedEquation('ude numerical', solve_buffer, solve_buffer_deriv, [c23, c24], params={'volume_tank': 1200, 't_tank_initial': t_buffer})
    nw.add_ude(ude)
    nw.set_attr(m_range=[.1, 100])
    nw.solve('design')

    return nw

q_owd = 9.37
q_tud = 6.07
t_supply = 85
t_return_owd = 46
q_condenser = 10
t_buffer = 80.64

res5 = option_1_mode5_test(q_owd, q_tud, t_supply, t_return_tud, t_return_owd, q_condenser, t_buffer)
res5.print_results()

Answered by fwitte

Apr 27, 2024

Hi Paula,

thank you, that is very helpful. Simple issue actually:

calc_T (all the calc_... methods of Connections actually) will provide you SI values... So you are mixing °C with K here, which is why, the solver cannot find a solution (see the warning message on the remaining residual value). Go with something like this:

def simple_system_ude():

    sp= Splitter('Splitter')
    me = Merge('Merge')
    va = Valve('Valve')
    boiler = SimpleHeatExchanger('Boiler')
    buffer = SimpleHeatExchanger('Buffer')
    si = Sink('sink')
    so = Source('source')

    c0 = Connection(so, 'out1', sp, 'in1', label='0')
    c1 = Connection(sp, 'out1', va, 'in1', label='1')
    c1a = Connection(va, 'o…

View full answer

fwitte · 2024-04-19T14:12:37Z

fwitte
Apr 19, 2024
Maintainer

Hi and thanks for reaching out!

Can you post a minimal working example where you want to replace some condition with the ude (an example, where the simulation works without the ude inserted and tell me which condition should be replaced by the ude). That would help a lot trouble shooting the issue.

Thanks a lot, have a nice weekend

Francesco

6 replies

fwitte Apr 24, 2024
Maintainer

Hi thanks for the answer,

that is quite a lot to process, I'll have to find some spare time for that. Or, can you send me something very small (not more than 5 connections?)? Or we could have a short online meeting, maybe next Monday?

Best

Francesco

gregorlein Apr 26, 2024
Author

Hi Francesco,
I tried to break it down into a simple system. I have three equations with three unknowns (marked in red). Equation I and II are solved within tespy and for the third equation, that is the analytical solution for a mixing temperature within a tank, I use an extra function 'solve_buffer_discharging'. My question is, can I incorporate equation III directly in tespy in order to reduce the computational time... If it's too complex to discuss here, a meeting on monday would work for me! Best, Paula

### set up the network
def simple_system(t_buffer):

    sp= Splitter('Splitter')
    me = Merge('Merge')
    va = Valve('Valve')
    boiler = SimpleHeatExchanger('Boiler')
    buffer = SimpleHeatExchanger('Buffer')
    si = Sink('sink')
    so = Source('source')

    c0 = Connection(so, 'out1', sp, 'in1', label='0')
    c1 = Connection(sp, 'out1', va, 'in1', label='1')
    c1a = Connection(va, 'out1', boiler, 'in1', label = '1a')
    c2 = Connection(boiler, 'out1', me, 'in1', label='2')
    c3 = Connection(buffer, 'out1', me, 'in2', label='3')
    c4 = Connection(sp, 'out2', buffer, 'in1', label='4')
    c5 = Connection(me, 'out1', si, 'in1', label='5')

    nw = Network(T_unit='C', p_unit='bar', h_unit='kJ / kg', m_unit = 'kg / s', v_unit = 'm3 / h')
    nw.add_conns(c0, c1, c1a, c2, c3, c4, c5)

    boiler.set_attr(pr=0.98)
    buffer.set_attr(pr=0.98)

    c0.set_attr(fluid={'water': 1}, T=50, p=1, v=300)
    c2.set_attr(T=95)
    c5.set_attr(T=90)
    c3.set_attr(T=t_buffer)
    nw.set_attr(iterinfo=False)
    nw.solve('design')

    return nw

t_init = 80 # initial temperature in buffer
tank_volume = 1200 # m³
delta_t = 10 # guess by how much buffer will cool down

### define iteration function to solve mixing temperature in buffer as well as tespy model
def solve_buffer_discharging(t_into_buffer, tol=1e-6, max_iter=10):
    k1 = t_init - t_into_buffer
    t_out_guess = t_init - delta_t # cool down buffer

    for _ in range(max_iter):
        result = simple_system(t_out_guess)         # solve system with temperature guess
        flow = abs(result.results['Connection'].loc['4'].v)  # result for flow with guessed temperature
        t_out = t_into_buffer + k1 * np.exp(-flow / tank_volume) # analytical solution for mixing temperature in buffer

        # Check convergence
        diff = abs(t_out_guess - t_out)
        if diff < tol:
            return result, t_out
        
        # if no convergence update guess
        t_out_guess = t_out
    
    print('no solution found')

### call function
t_into_buffer = 50 # temperature into buffer is temperatur in c0
result, t_out = solve_buffer_discharging(t_into_buffer)
result.print_results()

gregorlein Apr 26, 2024
Author

as an update, here is my idea for the network with an UDE. The code is working, but the result for the temperature in connection 3 doesn't seem correct to me...

def simple_system_ude():

    sp= Splitter('Splitter')
    me = Merge('Merge')
    va = Valve('Valve')
    boiler = SimpleHeatExchanger('Boiler')
    buffer = SimpleHeatExchanger('Buffer')
    si = Sink('sink')
    so = Source('source')

    c0 = Connection(so, 'out1', sp, 'in1', label='0')
    c1 = Connection(sp, 'out1', va, 'in1', label='1')
    c1a = Connection(va, 'out1', boiler, 'in1', label = '1a')
    c2 = Connection(boiler, 'out1', me, 'in1', label='2')
    c3 = Connection(buffer, 'out1', me, 'in2', label='3')
    c4 = Connection(sp, 'out2', buffer, 'in1', label='4')
    c5 = Connection(me, 'out1', si, 'in1', label='5')

    nw = Network(T_unit='C', p_unit='bar', h_unit='kJ / kg', m_unit = 'kg / s', v_unit = 'm3 / h')
    nw.add_conns(c0, c1, c1a, c2, c3, c4, c5)

    boiler.set_attr(pr=0.98)
    buffer.set_attr(pr=0.98)

    c0.set_attr(fluid={'water': 1}, T=50, p=1, v=300)
    c2.set_attr(T=95)
    c5.set_attr(T=90)
    #c3.set_attr(T=t_buffer)

    def my_ude(ude):
        t_init = ude.params['t_init']
        tank_volume = ude.params['tank_volume']
        v3 = (ude.conns[0].m.val_SI/1000)*3600
        return (ude.conns[1].calc_T() + (t_init - ude.conns[1].calc_T())*np.exp(-v3/tank_volume)-ude.conns[0].calc_T())
    
    def my_ude_deriv(ude):
        c3 = ude.conns[0]
        c4 = ude.conns[1]
        if c4.h.is_var:
            ude.jacobian[c4.h.J_col] = ude.numeric_deriv('h', c4)
        if c3.h.is_var:
            ude.jacobian[c3.h.J_col] = ude.numeric_deriv('h', c3)

    ude = UserDefinedEquation('my ude', my_ude, my_ude_deriv, [c3, c4], params={'t_init': 80, 'tank_volume': 1200})
    nw.add_ude(ude)
    nw.set_attr(m_range=[.1, 100])  # stabilize algorithm
    nw.set_attr(iterinfo=False)
    nw.solve('design')

    return nw

fwitte Apr 27, 2024
Maintainer

Hi Paula,

thank you, that is very helpful. Simple issue actually:

calc_T (all the calc_... methods of Connections actually) will provide you SI values... So you are mixing °C with K here, which is why, the solver cannot find a solution (see the warning message on the remaining residual value). Go with something like this:

def simple_system_ude():

    sp= Splitter('Splitter')
    me = Merge('Merge')
    va = Valve('Valve')
    boiler = SimpleHeatExchanger('Boiler')
    buffer = SimpleHeatExchanger('Buffer')
    si = Sink('sink')
    so = Source('source')

    c0 = Connection(so, 'out1', sp, 'in1', label='0')
    c1 = Connection(sp, 'out1', va, 'in1', label='1')
    c1a = Connection(va, 'out1', boiler, 'in1', label = '1a')
    c2 = Connection(boiler, 'out1', me, 'in1', label='2')
    c3 = Connection(buffer, 'out1', me, 'in2', label='3')
    c4 = Connection(sp, 'out2', buffer, 'in1', label='4')
    c5 = Connection(me, 'out1', si, 'in1', label='5')

    nw = Network(T_unit='C', p_unit='bar', h_unit='kJ / kg', m_unit = 'kg / s', v_unit = 'm3 / h')
    nw.add_conns(c0, c1, c1a, c2, c3, c4, c5)

    boiler.set_attr(pr=0.98)
    buffer.set_attr(pr=0.98)

    c0.set_attr(fluid={'water': 1}, T=50, p=1, v=300)
    c2.set_attr(T=95)
    c5.set_attr(T=90)
    # c3.set_attr(T=77.88)

    def my_ude(ude):
        c3 = ude.conns[0]
        c4 = ude.conns[1]
        t_init = ude.params['t_init']
        tank_volume = ude.params['tank_volume']
        v3 = c3.m.val_SI * c3.calc_vol() * 3600  # m3/h
        residual = (c3.calc_T() - 273.15) - ((c4.calc_T() - 273.15) + (t_init - (c4.calc_T() - 273.15)) * np.exp(-v3 / tank_volume))
        return residual

    def my_ude_deriv(ude):
        c3 = ude.conns[0]
        c4 = ude.conns[1]
        if c3.m.is_var:
            ude.jacobian[c3.m.J_col] = ude.numeric_deriv('m', c3)
        if c4.p.is_var:
            ude.jacobian[c4.p.J_col] = ude.numeric_deriv('p', c4)
        if c3.p.is_var:
            ude.jacobian[c3.p.J_col] = ude.numeric_deriv('p', c3)
        if c4.h.is_var:
            ude.jacobian[c4.h.J_col] = ude.numeric_deriv('h', c4)
        if c3.h.is_var:
            ude.jacobian[c3.h.J_col] = ude.numeric_deriv('h', c3)

    ude = UserDefinedEquation('my ude', my_ude, my_ude_deriv, [c3, c4], params={'t_init': 80, 'tank_volume': 1200})
    nw.add_ude(ude)
    nw.set_attr(m_range=[1, 100])  # stabilize algorithm
    nw.set_attr()
    nw.solve('design')
    nw.print_results()

    return nw

simple_system_ude()

I am getting nearly the same temperature as in your iterative implementation (diff of 0.02 °C).

Best

Francesco

Answer selected by gregorlein

gregorlein Apr 30, 2024
Author

Thanks for your help!! :)

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

UDE to solve buffer #505

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 1 comment 6 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

Select a reply

UDE to solve buffer #505

gregorlein Apr 18, 2024

Replies: 1 comment · 6 replies

fwitte Apr 19, 2024 Maintainer

fwitte Apr 24, 2024 Maintainer

gregorlein Apr 26, 2024 Author

gregorlein Apr 26, 2024 Author

fwitte Apr 27, 2024 Maintainer

gregorlein Apr 30, 2024 Author

gregorlein
Apr 18, 2024

Replies: 1 comment 6 replies

fwitte
Apr 19, 2024
Maintainer

fwitte Apr 24, 2024
Maintainer

gregorlein Apr 26, 2024
Author

gregorlein Apr 26, 2024
Author

fwitte Apr 27, 2024
Maintainer

gregorlein Apr 30, 2024
Author