optimizationminimizegekko

GEKKO Optimization of a DAE


I am trying to solve a DAE system using GEKKO. In the first part, I successfully solved the DAE. In the second part, I want to minimize one of my variables ('R'). Unfortunately, I cant find any solutions. Do you have any ideas what I did wrong? Here is my code:

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt


Tf = 150

T10 = 150.0
T20 = 150.0
T30 = 150.0
T40 = 150.0
T50 = 150.0
T60 = 150.0

Ta10 = Tf
Ta20 = Tf
Ta30 = Tf
Ta40 = Tf
Ta50 = Tf
Ta60 = Tf
T_f0 = Tf+1
Tw20 = 61

# Parameters:
# Stahl
m_steel = 7500 / 3600  # kg
cp = 500  # J/kg/K
# Luft
m_air = 100000 / 3600
cp_air = 10  # J/kg/K
alpha = 4500  # J/s
# Wasser
m_w = 0.15      # Strom durch HX


cp_w = 4200     # J/kg/K
# HX
UA = 120    # Wärmeübergangskoeffizient
Tw00 = 60   # Zulauftempratur


m = GEKKO()     # create GEKKO model; "remote=False" for connection issues
R = m.MV(0.0, lb=0.0, ub=0.9) # Reflux
# m_w = m.MV(0.02, lb = 0.02, ub=1)
t = np.linspace(0, 30, 100)
m.time = t

# Charge
Q_c = np.zeros(100)
Q_c[2:40] = 300000
Q_charge = m.Param(value=Q_c)

# Discharge
Q_d = np.zeros(100)
Q_d[60:80] = -1
Q_discharge = m.Param(value=Q_d)

T1, T2, T3, T4, T5, T6, Ta1, Ta2, Ta3, Ta4, Ta5, Ta6, T_f, Tw2, Tw1, Tw12, Q_Dis = m.Array(m.Var, 17)
T1.value = T10
T2.value = T20
T3.value = T30
T4.value = T40
T5.value = T50
T6.value = T60
Ta1.value = Ta10
Ta2.value = Ta20
Ta3.value = Ta30
Ta4.value = Ta40
Ta5.value = Ta50
Ta6.value = Ta60
T_f.value = T_f0
Tw1.value = Tw20
Tw2.value = Tw20
Tw12.value = Tw20
Q_Dis.value = 0

m.Equations([   # Storage
             T1.dt() == (alpha * (T_f - T1)) / (m_steel * cp),
             T2.dt() == (alpha * (Ta1 - T2) + Q_charge) / (m_steel * cp),
             T3.dt() == (alpha * (Ta2 - T3)) / (m_steel * cp),
             T4.dt() == (alpha * (Ta3 - T4)) / (m_steel * cp),
             T5.dt() == (alpha * (Ta4 - T5)) / (m_steel * cp),
             T6.dt() == (alpha * (Ta5 - T6)) / (m_steel * cp),
                # Calculate Temperature of air
             Ta1.dt() == -(alpha * (Ta1 - T1)) / (m_air * cp_air),
             Ta2.dt() == -(alpha * (Ta2 - T2)) / (m_air * cp_air),
             Ta3.dt() == -(alpha * (Ta3 - T3)) / (m_air * cp_air),
             Ta4.dt() == -(alpha * (Ta4 - T4)) / (m_air * cp_air),
             Ta5.dt() == -(alpha * (Ta5 - T5)) / (m_air * cp_air),
             Ta6.dt() == -(alpha * (Ta6 - T6)) / (m_air * cp_air),
                # HX
             T_f == Q_Dis / (m_air * cp_air) + Ta6,
             Tw2 == -Q_Dis / (m_w * cp_w) + Tw1,
             # Q_Dis == UA * (((T_f-Tw2) - (T_f - Tw1))/m.log((T_f-Tw2)/(T_f - Tw1))) * Q_discharge,
             Q_Dis == UA * (T_f - Tw1)/2 * Q_discharge,
                # Nodes of the water system
             Tw1 == ((m_w * (1 - (1-R))) * Tw12 + m_w * (1 - R) * Tw00) / ((m_w * (1 - (1 - R))) + m_w * (1 - R)),
             Tw12 == Tw2
             ])


# initialize with simulation
m.options.IMODE = 7
m.options.NODES = 3
m.solve()

# plot the prediction
plt.figure(figsize=(8, 5))

plt.subplot(3, 1, 1)
plt.plot(m.time, T1.value, 'r-', linewidth=3)
plt.plot(m.time, T2.value, 'b-', linewidth=3)
plt.plot(m.time, T3.value, 'g-', linewidth=3)
plt.plot(m.time, T4.value, 'r--', linewidth=3)
plt.plot(m.time, T5.value, 'b--', linewidth=3)
plt.plot(m.time, T6.value, 'g--', linewidth=3)
plt.plot(m.time, Tw2.value, 'b--', linewidth=4)
plt.plot(m.time, T_f.value, 'g-', linewidth=4)
plt.plot(m.time, Tw1.value, 'b:', linewidth=1)

plt.subplot(3, 1, 2)
plt.plot(m.time, Ta1.value, 'r-', linewidth=3)
plt.plot(m.time, Ta2.value, 'b-', linewidth=3)
plt.plot(m.time, Ta3.value, 'g-', linewidth=3)
plt.plot(m.time, Ta4.value, 'r--', linewidth=3)
plt.plot(m.time, Ta5.value, 'b--', linewidth=3)
plt.plot(m.time, Ta6.value, 'g--', linewidth=3)

plt.subplot(3, 1, 3)
plt.plot(m.time, Q_Dis.value, 'r-', linewidth=3)

plt.show()

# Optimize
m.options.IMODE = 6
# Tw2.UPPER = 250
# Tw2.LOWER = 110
R.Status = 1
m.options.SOLVER = 3
m.options.TIME_SHIFT = 0
T1.value = T1.value.value
T2.value = T2.value.value
T3.value = T3.value.value
T4.value = T4.value.value
T5.value = T5.value.value
T6.value = T6.value.value
Ta1.value = Ta1.value.value
Ta2.value = Ta2.value.value
Ta3.value = Ta3.value.value
Ta4.value = Ta4.value.value
Ta5.value = Ta5.value.value
Ta6.value = Ta6.value.value
T_f.value = T_f.value.value
Tw1.value = Tw1.value.value
Tw2.value = Tw2.value.value
Tw12.value = Tw12.value.value
Q_Dis.value = Q_Dis.value.value

m.Minimize(R)
m.solve(disp=True)

plt.figure(figsize=(8, 5))

plt.subplot(3, 1, 1)
plt.plot(m.time, T1.value, 'r-', linewidth=3)
plt.plot(m.time, T2.value, 'b-', linewidth=3)
plt.plot(m.time, T3.value, 'g-', linewidth=3)
plt.plot(m.time, T4.value, 'r--', linewidth=3)
plt.plot(m.time, T5.value, 'b--', linewidth=3)
plt.plot(m.time, T6.value, 'g--', linewidth=3)
plt.plot(m.time, Tw2.value, 'b--', linewidth=4)
plt.plot(m.time, T_f.value, 'g-', linewidth=4)
plt.plot(m.time, Tw1.value, 'b:', linewidth=1)

plt.subplot(3, 1, 2)
plt.plot(m.time, Ta1.value, 'r-', linewidth=3)
plt.plot(m.time, Ta2.value, 'b-', linewidth=3)
plt.plot(m.time, Ta3.value, 'g-', linewidth=3)
plt.plot(m.time, Ta4.value, 'r--', linewidth=3)
plt.plot(m.time, Ta5.value, 'b--', linewidth=3)
plt.plot(m.time, Ta6.value, 'g--', linewidth=3)

plt.subplot(3, 1, 3)
plt.plot(m.time, R.value, 'r-', linewidth=3)

plt.show()

Kind Regars Cassian

I tried to simplifie my code as you can see here:

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt


Tf = 150

T10 = 150.0
T20 = 150.0
T30 = 150.0
T40 = 150.0
T50 = 150.0
T60 = 150.0

Ta10 = Tf
Ta20 = Tf
Ta30 = Tf
Ta40 = Tf
Ta50 = Tf
Ta60 = Tf
T_f0 = Tf+1
Tw20 = 61

# Parameters:
# Stahl
m_steel = 7500 / 3600  # kg
cp = 500  # J/kg/K
# Luft
m_air = 100000 / 3600
cp_air = 10  # J/kg/K
alpha = 4500  # J/s
# Wasser
# R = 0.3         # Reflux
m_w = 0.15      # Strom durch HX
# m_w00 = 0.15    # Zustrom
# m_w22 = m_w * R     # Abstrom
# m_w12 = m_w * (1 - R)     # Reflux Strom

cp_w = 4200     # J/kg/K
# HX
UA = 120    # Wärmeübergangskoeffizient
Tw00 = 60   # Zulauftempratur


m = GEKKO()     # create GEKKO model; "remote=False" for connection issues
R = m.MV(0.5, lb=0.49, ub=0.6)
# t = np.linspace(0, 30, 100)
t = np.linspace(0, 30, 50)
m.time = t

# Charge
# Q_c = np.zeros(100)
# Q_c[2:40] = 300000

Q_c = np.zeros(50)
Q_c[1:20] = 300000
Q_charge = m.Param(value=Q_c)

# Discharge
# Q_d = np.zeros(100)
# Q_d[60:80] = -1
Q_d = np.zeros(50)
Q_d[30:40] = -1
Q_discharge = m.Param(value=Q_d)

# T1, T2, T3, T4, T5, T6, Ta1, Ta2, Ta3, Ta4, Ta5, Ta6, T_f, Tw2, Tw1, Tw12, Q_Dis = m.Array(m.Var, 17)
T1, T6, Ta1, Ta6, T_f, Tw2, Tw1, Q_Dis = m.Array(m.Var, 8)
T1.value = T10
# T2.value = T20
# T3.value = T30
# T4.value = T40
# T5.value = T50
T6.value = T60
Ta1.value = Ta10
# Ta2.value = Ta20
# Ta3.value = Ta30
# Ta4.value = Ta40
# Ta5.value = Ta50
Ta6.value = Ta60
T_f.value = T_f0
Tw1.value = Tw20
Tw2.value = Tw20
# Tw12.value = Tw20
Q_Dis.value = 0

m.Equations([   # Storage
             T_f / 1000 == (Q_Dis / (m_air * cp_air) + Ta6) / 1000,
             Tw2 / 1000 == (-Q_Dis / (m_w * cp_w) + Tw1) / 1000,
             # Q_Dis == UA * (((T_f-Tw2) - (T_f - Tw1))/m.log((T_f-Tw2)/(T_f - Tw1))) * Q_discharge,
             Q_Dis / 10000 ==  UA * (T_f - Tw1) / 2 * Q_discharge / 10000,
                # Nodes of the water system
             Tw1 / 1000 == ((m_w * (R)) * Tw2 + m_w * (1 - R) * Tw00) / ((m_w * (R)) + m_w * (1 - R)) / 1000,
             # Tw1 == R * Tw2 + R * Tw00 + Tw00,
             # Tw12 / 1000 == Tw2 / 1000,
             T1.dt() / 1000 == (alpha * (T_f - T1)) / (m_steel * cp) / 1000,
             # T2.dt() == (alpha * (Ta1 - T2) + Q_charge) / (m_steel * cp),
             # T3.dt() == (alpha * (Ta2 - T3)) / (m_steel * cp),
             # T4.dt() == (alpha * (Ta3 - T4)) / (m_steel * cp),
             # T5.dt() == (alpha * (Ta4 - T5)) / (m_steel * cp),
             T6.dt() / 1000 == (alpha * (Ta1 - T6)) / (m_steel * cp) / 1000,
                # Calculate Temperature of air
             Ta1.dt() / 1000 == -(alpha * (Ta1 - T1)) / (m_air * cp_air) / 1000,
             # Ta2.dt() == -(alpha * (Ta2 - T2)) / (m_air * cp_air),
             # Ta3.dt() == -(alpha * (Ta3 - T3)) / (m_air * cp_air),
             # Ta4.dt() == -(alpha * (Ta4 - T4)) / (m_air * cp_air),
             # Ta5.dt() == -(alpha * (Ta5 - T5)) / (m_air * cp_air),
             Ta6.dt() / 1000 == -(alpha * (Ta6 - T6)) / (m_air * cp_air) / 1000
             ])


# initialize with simulation
m.options.IMODE = 7
m.options.NODES = 3
m.solve()

# plot the prediction
plt.figure(figsize=(8, 5))

plt.subplot(3, 1, 1)
plt.plot(m.time, T1.value, 'r-', linewidth=3)
# plt.plot(m.time, T2.value, 'b-', linewidth=3)
# plt.plot(m.time, T3.value, 'g-', linewidth=3)
# plt.plot(m.time, T4.value, 'r--', linewidth=3)
# plt.plot(m.time, T5.value, 'b--', linewidth=3)
plt.plot(m.time, T6.value, 'g--', linewidth=3)
plt.plot(m.time, Tw2.value, 'b--', linewidth=4)
plt.plot(m.time, T_f.value, 'g-', linewidth=4)
plt.plot(m.time, Tw1.value, 'b:', linewidth=1)

plt.subplot(3, 1, 2)
plt.plot(m.time, Ta1.value, 'r-', linewidth=3)
# plt.plot(m.time, Ta2.value, 'b-', linewidth=3)
# plt.plot(m.time, Ta3.value, 'g-', linewidth=3)
# plt.plot(m.time, Ta4.value, 'r--', linewidth=3)
# plt.plot(m.time, Ta5.value, 'b--', linewidth=3)
plt.plot(m.time, Ta6.value, 'g--', linewidth=3)

plt.subplot(3, 1, 3)
plt.plot(m.time, Q_Dis.value, 'r-', linewidth=3)

plt.show()

# Optimize
m.options.IMODE = 6
Q_d = np.zeros(50)
Q_d[30:40] = -1
Q_discharge = m.Param(value=Q_d)

Q_c = np.zeros(50)
Q_c[1:20] = 300000
Q_charge = m.Param(value=Q_c)
# Tw1.UPPER = 1250
# Tw2.LOWER = 10
R.Status = 1
m.options.SOLVER = 3
m.options.TIME_SHIFT = 0
T1.value = T1.value.value
# T2.value = T2.value.value
# T3.value = T3.value.value
# T4.value = T4.value.value
# T5.value = T5.value.value
T6.value = T6.value.value
Ta1.value = Ta1.value.value
# Ta2.value = Ta2.value.value
# Ta3.value = Ta3.value.value
# Ta4.value = Ta4.value.value
# Ta5.value = Ta5.value.value
Ta6.value = Ta6.value.value
T_f.value = T_f.value.value
Tw1.value = Tw1.value.value
Tw2.value = Tw2.value.value
# Tw12.value = Tw12.value.value
Q_Dis.value = 0 # Q_Dis.value.value

m.Minimize(R)
m.solve(disp=True)

plt.figure(figsize=(8, 5))

plt.subplot(3, 1, 1)
plt.plot(m.time, T1.value, 'r-', linewidth=3)
# plt.plot(m.time, T2.value, 'b-', linewidth=3)
# plt.plot(m.time, T3.value, 'g-', linewidth=3)
# plt.plot(m.time, T4.value, 'r--', linewidth=3)
# plt.plot(m.time, T5.value, 'b--', linewidth=3)
plt.plot(m.time, T6.value, 'g--', linewidth=3)
plt.plot(m.time, Tw2.value, 'b--', linewidth=4)
plt.plot(m.time, T_f.value, 'g-', linewidth=4)
plt.plot(m.time, Tw1.value, 'b:', linewidth=1)

plt.subplot(3, 1, 2)
plt.plot(m.time, Ta1.value, 'r-', linewidth=3)
# plt.plot(m.time, Ta2.value, 'b-', linewidth=3)
# plt.plot(m.time, Ta3.value, 'g-', linewidth=3)
# plt.plot(m.time, Ta4.value, 'r--', linewidth=3)
# plt.plot(m.time, Ta5.value, 'b--', linewidth=3)
plt.plot(m.time, Ta6.value, 'g--', linewidth=3)

plt.subplot(3, 1, 3)
plt.plot(m.time, R.value, 'r-', linewidth=3)

plt.show()

This code produces a solution, but not the one I expect or want haha. Some of my variables (Tw1 and Tw2) go against -∞ and the rest against +∞.


Solution

  • The simulation solves successfully, but the result is an unstable system.

    unstable

    The only objective is the minimize R. Try adding another objective such as:

    m.Minimize((Tw1-5)**2)
    

    and also widen the bounds for R such as R = m.MV(0.5, lb=0.2, ub=1.0). This gives a stable solution and an objective to guide the selection of R.

    stable solution

    import numpy as np
    from gekko import GEKKO
    import matplotlib.pyplot as plt
    
    Tf = 150
    
    T10 = 150.0
    T20 = 150.0
    T30 = 150.0
    T40 = 150.0
    T50 = 150.0
    T60 = 150.0
    
    Ta10 = Tf
    Ta20 = Tf
    Ta30 = Tf
    Ta40 = Tf
    Ta50 = Tf
    Ta60 = Tf
    T_f0 = Tf+1
    Tw20 = 61
    
    # Parameters:
    # Stahl
    m_steel = 7500 / 3600  # kg
    cp = 500  # J/kg/K
    # Luft
    m_air = 100000 / 3600
    cp_air = 10  # J/kg/K
    alpha = 4500  # J/s
    # Wasser
    # R = 0.3         # Reflux
    m_w = 0.15      # Strom durch HX
    # m_w00 = 0.15    # Zustrom
    # m_w22 = m_w * R     # Abstrom
    # m_w12 = m_w * (1 - R)     # Reflux Strom
    
    cp_w = 4200     # J/kg/K
    # HX
    UA = 120    # Wärmeübergangskoeffizient
    Tw00 = 60   # Zulauftempratur
    
    
    m = GEKKO()     # create GEKKO model; "remote=False" for connection issues
    R = m.MV(0.5, lb=0.2, ub=1.0)
    # t = np.linspace(0, 30, 100)
    t = np.linspace(0, 30, 50)
    m.time = t
    
    # Charge
    # Q_c = np.zeros(100)
    # Q_c[2:40] = 300000
    
    Q_c = np.zeros(50)
    Q_c[1:20] = 300000
    Q_charge = m.Param(value=Q_c)
    
    # Discharge
    # Q_d = np.zeros(100)
    # Q_d[60:80] = -1
    Q_d = np.zeros(50)
    Q_d[30:40] = -1
    Q_discharge = m.Param(value=Q_d)
    
    # T1, T2, T3, T4, T5, T6, Ta1, Ta2, Ta3, Ta4, Ta5, Ta6, T_f, Tw2, Tw1, Tw12, Q_Dis = m.Array(m.Var, 17)
    T1, T6, Ta1, Ta6, T_f, Tw2, Tw1, Q_Dis = m.Array(m.Var, 8)
    T1.value = T10
    # T2.value = T20
    # T3.value = T30
    # T4.value = T40
    # T5.value = T50
    T6.value = T60
    Ta1.value = Ta10
    # Ta2.value = Ta20
    # Ta3.value = Ta30
    # Ta4.value = Ta40
    # Ta5.value = Ta50
    Ta6.value = Ta60
    T_f.value = T_f0
    Tw1.value = Tw20
    Tw2.value = Tw20
    # Tw12.value = Tw20
    Q_Dis.value = 0
    
    m.Equations([   # Storage
                 T_f / 1000 == (Q_Dis / (m_air * cp_air) + Ta6) / 1000,
                 Tw2 / 1000 == (-Q_Dis / (m_w * cp_w) + Tw1) / 1000,
                 # Q_Dis == UA * (((T_f-Tw2) - (T_f - Tw1))/m.log((T_f-Tw2)/(T_f - Tw1))) * Q_discharge,
                 Q_Dis / 10000 ==  UA * (T_f - Tw1) / 2 * Q_discharge / 10000,
                    # Nodes of the water system
                 Tw1 / 1000 == ((m_w * (R)) * Tw2 + m_w * (1 - R) * Tw00) / ((m_w * (R)) + m_w * (1 - R)) / 1000,
                 # Tw1 == R * Tw2 + R * Tw00 + Tw00,
                 # Tw12 / 1000 == Tw2 / 1000,
                 T1.dt() / 1000 == (alpha * (T_f - T1)) / (m_steel * cp) / 1000,
                 # T2.dt() == (alpha * (Ta1 - T2) + Q_charge) / (m_steel * cp),
                 # T3.dt() == (alpha * (Ta2 - T3)) / (m_steel * cp),
                 # T4.dt() == (alpha * (Ta3 - T4)) / (m_steel * cp),
                 # T5.dt() == (alpha * (Ta4 - T5)) / (m_steel * cp),
                 T6.dt() / 1000 == (alpha * (Ta1 - T6)) / (m_steel * cp) / 1000,
                    # Calculate Temperature of air
                 Ta1.dt() / 1000 == -(alpha * (Ta1 - T1)) / (m_air * cp_air) / 1000,
                 # Ta2.dt() == -(alpha * (Ta2 - T2)) / (m_air * cp_air),
                 # Ta3.dt() == -(alpha * (Ta3 - T3)) / (m_air * cp_air),
                 # Ta4.dt() == -(alpha * (Ta4 - T4)) / (m_air * cp_air),
                 # Ta5.dt() == -(alpha * (Ta5 - T5)) / (m_air * cp_air),
                 Ta6.dt() / 1000 == -(alpha * (Ta6 - T6)) / (m_air * cp_air) / 1000
                 ])
    
    
    # initialize with simulation
    m.options.IMODE = 7
    m.options.NODES = 3
    m.solve()
    
    # plot the prediction
    plt.figure(figsize=(8, 5))
    
    plt.subplot(3, 1, 1)
    plt.plot(m.time, T1.value, 'r-', linewidth=3)
    # plt.plot(m.time, T2.value, 'b-', linewidth=3)
    # plt.plot(m.time, T3.value, 'g-', linewidth=3)
    # plt.plot(m.time, T4.value, 'r--', linewidth=3)
    # plt.plot(m.time, T5.value, 'b--', linewidth=3)
    plt.plot(m.time, T6.value, 'g--', linewidth=3)
    plt.plot(m.time, Tw2.value, 'b--', linewidth=4)
    plt.plot(m.time, T_f.value, 'g-', linewidth=4)
    plt.plot(m.time, Tw1.value, 'b:', linewidth=1)
    
    plt.subplot(3, 1, 2)
    plt.plot(m.time, Ta1.value, 'r-', linewidth=3)
    # plt.plot(m.time, Ta2.value, 'b-', linewidth=3)
    # plt.plot(m.time, Ta3.value, 'g-', linewidth=3)
    # plt.plot(m.time, Ta4.value, 'r--', linewidth=3)
    # plt.plot(m.time, Ta5.value, 'b--', linewidth=3)
    plt.plot(m.time, Ta6.value, 'g--', linewidth=3)
    
    plt.subplot(3, 1, 3)
    plt.plot(m.time, Q_Dis.value, 'r-', linewidth=3)
    
    plt.show()
    
    # Optimize
    m.options.IMODE = 6
    Q_d = np.zeros(50)
    Q_d[30:40] = -1
    Q_discharge = m.Param(value=Q_d)
    
    Q_c = np.zeros(50)
    Q_c[1:20] = 300000
    Q_charge = m.Param(value=Q_c)
    # Tw1.UPPER = 1250
    # Tw2.LOWER = 10
    R.Status = 1
    m.options.SOLVER = 3
    m.options.TIME_SHIFT = 0
    T1.value = T1.value.value
    # T2.value = T2.value.value
    # T3.value = T3.value.value
    # T4.value = T4.value.value
    # T5.value = T5.value.value
    T6.value = T6.value.value
    Ta1.value = Ta1.value.value
    # Ta2.value = Ta2.value.value
    # Ta3.value = Ta3.value.value
    # Ta4.value = Ta4.value.value
    # Ta5.value = Ta5.value.value
    Ta6.value = Ta6.value.value
    T_f.value = T_f.value.value
    Tw1.value = Tw1.value.value
    Tw2.value = Tw2.value.value
    # Tw12.value = Tw12.value.value
    Q_Dis.value = 0 # Q_Dis.value.value
    
    m.Minimize(R)
    m.Minimize((Tw1-5)**2)
    m.solve(disp=True)
    
    plt.figure(figsize=(8, 5))
    
    plt.subplot(3, 1, 1)
    plt.plot(m.time, T1.value, 'r-', linewidth=3)
    # plt.plot(m.time, T2.value, 'b-', linewidth=3)
    # plt.plot(m.time, T3.value, 'g-', linewidth=3)
    # plt.plot(m.time, T4.value, 'r--', linewidth=3)
    # plt.plot(m.time, T5.value, 'b--', linewidth=3)
    plt.plot(m.time, T6.value, 'g--', linewidth=3)
    plt.plot(m.time, Tw2.value, 'b--', linewidth=4)
    plt.plot(m.time, T_f.value, 'g-', linewidth=4)
    plt.plot(m.time, Tw1.value, 'b:', linewidth=1)
    
    plt.subplot(3, 1, 2)
    plt.plot(m.time, Ta1.value, 'r-', linewidth=3)
    # plt.plot(m.time, Ta2.value, 'b-', linewidth=3)
    # plt.plot(m.time, Ta3.value, 'g-', linewidth=3)
    # plt.plot(m.time, Ta4.value, 'r--', linewidth=3)
    # plt.plot(m.time, Ta5.value, 'b--', linewidth=3)
    plt.plot(m.time, Ta6.value, 'g--', linewidth=3)
    
    plt.subplot(3, 1, 3)
    plt.plot(m.time, R.value, 'r-', linewidth=3)
    
    plt.show()