<reponame>dhill2522/ChE436-project<gh_stars>0
import tclab  # pip install tclab
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# FOPDT model
Kp = 0.5      # degC/%
tauP = 120.0  # seconds
thetaP = 10   # seconds (integer)
Tss = 23      # degC (ambient temperature)
Qss = 0       # % heater

def heat(x, t, Q):
    '''Non-linear heater energy balance'''
    # Parameters
    Ta = 23 + 273.15   # K
    U = 10.0           # W/m^2-K
    m = 4.0/1000.0     # kg
    Cp = 0.5 * 1000.0  # J/kg-K
    A = 12.0 / 100.0**2  # Area in m^2
    alpha = 0.01       # W / % heater
    eps = 0.9          # Emissivity
    sigma = 5.67e-8    # Stefan-Boltzman

    # Temperature State
    T = x[0]

    # Nonlinear Energy Balance
    dTdt = (1.0/(m*Cp))*(U*A*(Ta-T)
                         + eps * sigma * A * (Ta**4 - T**4)
                         + alpha*Q)
    return dTdt

def save_txt(t, u1, u2, y1, y2, sp1, sp2):
    '''save data to a text file'''
    data = np.vstack((t, u1, u2, y1, y2, sp1, sp2))  # vertical stack
    data = data.T                 # transpose data
    top = 'Time (sec), Heater 1 (%), Heater 2 (%), ' \
        + 'Temperature 1 (degC), Temperature 2 (degC), ' \
        + 'Set Point 1 (degC), Set Point 2 (degC)'
    np.savetxt('data.txt', data, delimiter=',', header=top, comments='')


# Connect to Arduino
a = tclab.TCLab()

# Turn LED on
print('LED On')
a.LED(100)

# Run time in minutes
run_time = 30.0

# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)

# Temperature (K)
Tsp1 = np.ones(loops) * 23.0  # set point (degC)
T1 = np.ones(loops) * a.T1  # measured T (degC)

Tsp2 = np.ones(loops) * 23.0  # set point (degC)
T2 = np.ones(loops) * a.T2  # measured T (degC)

# Predictions
Tp = np.ones(loops) * a.T1
error_eb = np.zeros(loops)
Tpl = np.ones(loops) * a.T1
error_fopdt = np.zeros(loops)

# impulse tests (0 - 100%)
Q1 = np.ones(loops) * 0.0
Q2 = np.ones(loops) * 0.0
Q1[10:900] = 30.0  # step up for 100 sec
Q1[900:1800] = 20.0  # step up for 100 sec

print('Running Main Loop. Ctrl-C to end.')
print('  Time   Q1     Q2    T1     T2')
print('{:6.1f} {:6.2f} {:6.2f} {:6.2f} {:6.2f}'.format(tm[0],
                                                       Q1[0],
                                                       Q2[0],
                                                       T1[0],
                                                       T2[0]))

# Create plot
plt.figure(figsize=(10, 7))
plt.ion()
plt.show()

# Main Loop
start_time = time.time()
prev_time = start_time
try:
    for i in range(1, loops):
        # Sleep time
        sleep_max = 1.0
        sleep = sleep_max - (time.time() - prev_time)
        if sleep >= 0.01:
            time.sleep(sleep-0.01)
        else:
            time.sleep(0.01)

        # Record time and change in time
        t = time.time()
        dt = t - prev_time
        prev_time = t
        tm[i] = t - start_time

        # Read temperatures in Kelvin
        T1[i] = a.T1
        T2[i] = a.T2

        # Simulate one time step with Energy Balance
        Tnext = odeint(heat, Tp[i-1]+273.15, [0, dt], args=(Q1[i-1],))
        Tp[i] = Tnext[1]-273.15
        error_eb[i] = error_eb[i-1] + abs(Tp[i]-T1[i])

        # Simulate one time step with FOPDT model
        z = np.exp(-dt/tauP)
        Tpl[i] = (Tpl[i-1]-Tss) * z \
            + (Q1[max(0, i-int(thetaP)-1)]-Qss)*(1-z)*Kp \
            + Tss
        error_fopdt[i] = error_fopdt[i-1] + abs(Tpl[i]-T1[i])

        # Write output (0-100)
        a.Q1(Q1[i])
        a.Q2(Q2[i])

        # Print line of data
        print('{:6.1f} {:6.2f} {:6.2f} {:6.2f} {:6.2f}'.format(tm[i],
                                                               Q1[i],
                                                               Q2[i],
                                                               T1[i],
                                                               T2[i]))

        # Plot
        plt.clf()
        ax = plt.subplot(3, 1, 1)
        ax.grid()
        plt.plot(tm[0:i], T1[0:i], 'ro', label=r'$T_1$ measured')
        plt.plot(tm[0:i], Tp[0:i], 'k-', label=r'$T_1$ energy balance')
        plt.plot(tm[0:i], Tpl[0:i], 'g:', label=r'$T_1$ FOPDT')
        plt.plot(tm[0:i], T2[0:i], 'bx', label=r'$T_2$ measured')
        plt.ylabel('Temperature (degC)')
        plt.legend(loc=2)
        ax = plt.subplot(3, 1, 2)
        ax.grid()
        plt.plot(tm[0:i], error_eb[0:i], 'k-', label='Energy Balance')
        plt.plot(tm[0:i], error_fopdt[0:i], 'g:', label='Linear')
        plt.ylabel('Cumulative Error')
        plt.legend(loc='best')
        ax = plt.subplot(3, 1, 3)
        ax.grid()
        plt.plot(tm[0:i], Q1[0:i], 'r-', label=r'$Q_1$')
        plt.plot(tm[0:i], Q2[0:i], 'b:', label=r'$Q_2$')
        plt.ylabel('Heaters')
        plt.xlabel('Time (sec)')
        plt.legend(loc='best')
        plt.draw()
        plt.pause(0.05)

    # Turn off heaters
    a.Q1(0)
    a.Q2(0)
    # Save text file
    save_txt(tm[0:i], Q1[0:i], Q2[0:i], T1[0:i], T2[0:i], Tsp1[0:i], Tsp2[0:i])
    # Save figure
    plt.savefig('test_Models.eps', format='eps', dpi=1000)

# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
    # Disconnect from Arduino
    a.Q1(0)
    a.Q2(0)
    print('Shutting down')
    a.close()
    save_txt(tm[0:i], Q1[0:i], Q2[0:i], T1[0:i], T2[0:i], Tsp1[0:i], Tsp2[0:i])
    plt.savefig('test_Models.eps', format='eps', dpi=1000)

# Make sure serial connection still closes when there's an error
except:
    # Disconnect from Arduino
    a.Q1(0)
    a.Q2(0)
    print('Error: Shutting down')
    a.close()
    save_txt(tm[0:i], Q1[0:i], Q2[0:i], T1[0:i], T2[0:i], Tsp1[0:i], Tsp2[0:i])
    plt.savefig('test_Models.eps', format='eps', dpi=1000)
    raise

import tclab  # pip install tclab
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# FOPDT model
Kp = 0.5      # degC/%
tauP = 120.0  # seconds
thetaP = 10   # seconds (integer)
Tss = 23      # degC (ambient temperature)
Qss = 0       # % heater

def heat(x, t, Q):
    '''Non-linear heater energy balance'''
    # Parameters
    Ta = 23 + 273.15   # K
    U = 10.0           # W/m^2-K
    m = 4.0/1000.0     # kg
    Cp = 0.5 * 1000.0  # J/kg-K
    A = 12.0 / 100.0**2  # Area in m^2
    alpha = 0.01       # W / % heater
    eps = 0.9          # Emissivity
    sigma = 5.67e-8    # Stefan-Boltzman

    # Temperature State
    T = x[0]

    # Nonlinear Energy Balance
    dTdt = (1.0/(m*Cp))*(U*A*(Ta-T)
                         + eps * sigma * A * (Ta**4 - T**4)
                         + alpha*Q)
    return dTdt

def save_txt(t, u1, u2, y1, y2, sp1, sp2):
    '''save data to a text file'''
    data = np.vstack((t, u1, u2, y1, y2, sp1, sp2))  # vertical stack
    data = data.T                 # transpose data
    top = 'Time (sec), Heater 1 (%), Heater 2 (%), ' \
        + 'Temperature 1 (degC), Temperature 2 (degC), ' \
        + 'Set Point 1 (degC), Set Point 2 (degC)'
    np.savetxt('data.txt', data, delimiter=',', header=top, comments='')


# Connect to Arduino
a = tclab.TCLab()

# Turn LED on
print('LED On')
a.LED(100)

# Run time in minutes
run_time = 30.0

# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)

# Temperature (K)
Tsp1 = np.ones(loops) * 23.0  # set point (degC)
T1 = np.ones(loops) * a.T1  # measured T (degC)

Tsp2 = np.ones(loops) * 23.0  # set point (degC)
T2 = np.ones(loops) * a.T2  # measured T (degC)

# Predictions
Tp = np.ones(loops) * a.T1
error_eb = np.zeros(loops)
Tpl = np.ones(loops) * a.T1
error_fopdt = np.zeros(loops)

# impulse tests (0 - 100%)
Q1 = np.ones(loops) * 0.0
Q2 = np.ones(loops) * 0.0
Q1[10:900] = 30.0  # step up for 100 sec
Q1[900:1800] = 20.0  # step up for 100 sec

print('Running Main Loop. Ctrl-C to end.')
print('  Time   Q1     Q2    T1     T2')
print('{:6.1f} {:6.2f} {:6.2f} {:6.2f} {:6.2f}'.format(tm[0],
                                                       Q1[0],
                                                       Q2[0],
                                                       T1[0],
                                                       T2[0]))

# Create plot
plt.figure(figsize=(10, 7))
plt.ion()
plt.show()

# Main Loop
start_time = time.time()
prev_time = start_time
try:
    for i in range(1, loops):
        # Sleep time
        sleep_max = 1.0
        sleep = sleep_max - (time.time() - prev_time)
        if sleep >= 0.01:
            time.sleep(sleep-0.01)
        else:
            time.sleep(0.01)

        # Record time and change in time
        t = time.time()
        dt = t - prev_time
        prev_time = t
        tm[i] = t - start_time

        # Read temperatures in Kelvin
        T1[i] = a.T1
        T2[i] = a.T2

        # Simulate one time step with Energy Balance
        Tnext = odeint(heat, Tp[i-1]+273.15, [0, dt], args=(Q1[i-1],))
        Tp[i] = Tnext[1]-273.15
        error_eb[i] = error_eb[i-1] + abs(Tp[i]-T1[i])

        # Simulate one time step with FOPDT model
        z = np.exp(-dt/tauP)
        Tpl[i] = (Tpl[i-1]-Tss) * z \
            + (Q1[max(0, i-int(thetaP)-1)]-Qss)*(1-z)*Kp \
            + Tss
        error_fopdt[i] = error_fopdt[i-1] + abs(Tpl[i]-T1[i])

        # Write output (0-100)
        a.Q1(Q1[i])
        a.Q2(Q2[i])

        # Print line of data
        print('{:6.1f} {:6.2f} {:6.2f} {:6.2f} {:6.2f}'.format(tm[i],
                                                               Q1[i],
                                                               Q2[i],
                                                               T1[i],
                                                               T2[i]))

        # Plot
        plt.clf()
        ax = plt.subplot(3, 1, 1)
        ax.grid()
        plt.plot(tm[0:i], T1[0:i], 'ro', label=r'$T_1$ measured')
        plt.plot(tm[0:i], Tp[0:i], 'k-', label=r'$T_1$ energy balance')
        plt.plot(tm[0:i], Tpl[0:i], 'g:', label=r'$T_1$ FOPDT')
        plt.plot(tm[0:i], T2[0:i], 'bx', label=r'$T_2$ measured')
        plt.ylabel('Temperature (degC)')
        plt.legend(loc=2)
        ax = plt.subplot(3, 1, 2)
        ax.grid()
        plt.plot(tm[0:i], error_eb[0:i], 'k-', label='Energy Balance')
        plt.plot(tm[0:i], error_fopdt[0:i], 'g:', label='Linear')
        plt.ylabel('Cumulative Error')
        plt.legend(loc='best')
        ax = plt.subplot(3, 1, 3)
        ax.grid()
        plt.plot(tm[0:i], Q1[0:i], 'r-', label=r'$Q_1$')
        plt.plot(tm[0:i], Q2[0:i], 'b:', label=r'$Q_2$')
        plt.ylabel('Heaters')
        plt.xlabel('Time (sec)')
        plt.legend(loc='best')
        plt.draw()
        plt.pause(0.05)

    # Turn off heaters
    a.Q1(0)
    a.Q2(0)
    # Save text file
    save_txt(tm[0:i], Q1[0:i], Q2[0:i], T1[0:i], T2[0:i], Tsp1[0:i], Tsp2[0:i])
    # Save figure
    plt.savefig('test_Models.eps', format='eps', dpi=1000)

# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
    # Disconnect from Arduino
    a.Q1(0)
    a.Q2(0)
    print('Shutting down')
    a.close()
    save_txt(tm[0:i], Q1[0:i], Q2[0:i], T1[0:i], T2[0:i], Tsp1[0:i], Tsp2[0:i])
    plt.savefig('test_Models.eps', format='eps', dpi=1000)

# Make sure serial connection still closes when there's an error
except:
    # Disconnect from Arduino
    a.Q1(0)
    a.Q2(0)
    print('Error: Shutting down')
    a.close()
    save_txt(tm[0:i], Q1[0:i], Q2[0:i], T1[0:i], T2[0:i], Tsp1[0:i], Tsp2[0:i])
    plt.savefig('test_Models.eps', format='eps', dpi=1000)
    raise
# pip install tclab # FOPDT model # degC/% # seconds # seconds (integer) # degC (ambient temperature) # % heater Non-linear heater energy balance # Parameters # K # W/m^2-K # kg # J/kg-K # Area in m^2 # W / % heater # Emissivity # Stefan-Boltzman # Temperature State # Nonlinear Energy Balance save data to a text file # vertical stack # transpose data # Connect to Arduino # Turn LED on # Run time in minutes # Number of cycles # Temperature (K) # set point (degC) # measured T (degC) # set point (degC) # measured T (degC) # Predictions # impulse tests (0 - 100%) # step up for 100 sec # step up for 100 sec # Create plot # Main Loop # Sleep time # Record time and change in time # Read temperatures in Kelvin # Simulate one time step with Energy Balance # Simulate one time step with FOPDT model # Write output (0-100) # Print line of data # Plot # Turn off heaters # Save text file # Save figure # Allow user to end loop with Ctrl-C # Disconnect from Arduino # Make sure serial connection still closes when there's an error # Disconnect from Arduino