4.3. Four State Model#

4.3.1. Multivariable Heater/Sensor System#

Model equations

\[\begin{split} \begin{align} C^H_p\frac{dT_{H,1}}{dt} & = U_a(T_{amb} - T_{H,1}) + U_b(T_{H,2}-T_{H,1}) + U_c(T_{S,1} - T_{H,1}) + P_1u_1\\ C^S_p\frac{dT_{S,1}}{dt} & = U_c(T_{H,1} - T_{S,1}) \\ C^H_p\frac{dT_{H,2}}{dt} & = U_a(T_{amb} - T_{H,2}) + U_b(T_{H,1}-T_{H,2}) + U_c(T_{S,2} - T_{H,2}) + P_2u_2\\ C^S_p\frac{dT_{S,2}}{dt} & = U_c(T_{H,2} - T_{S,2}) \end{align} \end{split}\]

4.3.2. Deviation variables#

\[\begin{split} \begin{align} \frac{dT_{H,1}'}{dt} & = -(\frac{U_a+U_b+U_c}{C^H_p})T_{H,1}' + \frac{U_b}{C^H_p}T_{H,2}' + \frac{U_c}{C^H_p}T_{S,1}' + \frac{P_1}{C^H_p}u_1\\ \frac{dT_{S,1}'}{dt} & = \frac{U_c}{C^S_p}(T_{H,1}' - T_{S,1}') \\ \frac{dT_{H,2}'}{dt} & = -(\frac{U_a+U_b+U_c}{C^H_p})T_{H,2}' + \frac{U_b}{C^H_p}T_{H,1}' + \frac{U_c}{C^H_p}T_{S,2}' + \frac{P_2}{C^H_p}u_2\\ \frac{dT_{S,2}'}{dt} & = \frac{U_c}{C^S_p}(T_{H,2}' - T_{S,2}') \end{align} \end{split}\]

4.3.3. State space#

\[\begin{split} \begin{align} \left[\begin{array}{c}\frac{dT_{H,1}'}{dt} \\ \frac{dT_{S,1}'}{dt} \\ \frac{dT_{H,2}'}{dt} \\ \frac{dT_{S,2}'}{dt}\end{array}\right] & = \left[\begin{array}{cccc} -(\frac{U_a+U_b+U_c}{C^H_p}) & \frac{U_c}{C^H_p} & \frac{U_b}{C^H_p} & 0 \\ \frac{U_c}{C^S_p} & -\frac{U_c}{C^S_p} & 0 & 0 \\ \frac{U_b}{C^H_p} & 0 & -(\frac{U_a+U_b+U_c}{C^H_p}) & \frac{U_c}{C^H_p} \\ 0 & 0 & \frac{U_c}{C^H_p} & -\frac{U_c}{C^H_p} \end{array}\right] \left[\begin{array}{c}T_{H,1}' \\ T_{S,1}' \\ T_{H,2}' \\ T_{S,2}'\end{array}\right] + \left[\begin{array}{cc}\frac{P_1}{C_p} & 0 \\ 0 & 0 \\ 0 & \frac{P_2}{C_p} \\ 0 & 0 \end{array}\right] \left[\begin{array}{c}u_1 \\ u_2\end{array}\right] \end{align} \end{split}\]
%matplotlib inline
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from ipywidgets import interact

data = pd.read_csv('Step_Test_Data.csv').set_index('Time')[1:]
t = data.index
T1 = data['T1'].values
T2 = data['T2'].values

# known parameter values
P1 = 4
u1 = 0.5   # steady state value of u1 (fraction of total power)
P2 = 2
u2 = 0.0
T_ambient = 21.5

def compare(Ua, Ub, Uc, Cp_H, Cp_S):

    def deriv(T,t):
        T_H1, T_S1, T_H2, T_S2 = T
        dT_H1 = -(Ua + Ub + Uc)*T_H1/Cp_H + Ub*T_H2/Cp_H + Uc*T_S1/Cp_H + P1*u1/Cp_H
        dT_S1 = Uc*T_H1/Cp_S - Uc*T_S1/Cp_S
        dT_H2 = -(Ua + Ub + Uc)*T_H2/Cp_H + Ub*T_H1/Cp_H + Uc*T_S2/Cp_H + P2*u2/Cp_H
        dT_S2 = Uc*T_H2/Cp_S - Uc*T_S2/Cp_S
        return [dT_H1, dT_S1, dT_H2, dT_S2]
    T = odeint(deriv, [0,0,0,0], t)

    plt.figure(figsize=(12,6))
    plt.subplot(2,1,1)
    plt.plot(t, T[:,1] + T_ambient, t, T[:,3] + T_ambient, t, T1, t, T2)
    plt.xlabel('Time / seconds')
    plt.ylabel('Temperature / °C')
    plt.grid()
    plt.subplot(2,1,2)
    plt.plot(t, T[:,1]-T1+T_ambient)
    plt.plot(t, T[:,3]-T2+T_ambient)
    plt.grid

# parameter values and units
P1 = 4                 # watts
P2 = 2                 # watts
Ua = 0.044             # watts/deg C
Ub = 0.018             # watts/deg C
Cp = 7                 # joules/deg C

w = interact(compare, Ua=(0,0.06,0.001), Ub=(0,0.06,0.001), Uc=(0,0.06,0.001), Cp_H=(3,11,0.01), Cp_S = (0.1,2,.01))
import numpy as np

Ua = 0.043
Ub = 0.022
Uc = 0.036
Cp_H = 6.38
Cp_S = 0.98

A = [[-(Ua+Ub+Uc)/Cp_H, Uc/Cp_H, Ub/Cp_H, 0],
     [Uc/Cp_S, -Uc/Cp_S, 0, 0],
     [Ub/Cp_H, 0, -(Ua+Ub+Uc)/Cp_H, Uc/Cp_H],
     [0, 0, Uc/Cp_H, -Uc/Cp_H]
    ]

evals, evecs = np.linalg.eig(A)

1/evals
array([ -22.64294396,  -52.59946787, -124.34799519, -354.43970971])
def compare(Ua, Ub, Uc, Cp_H, Cp_S):

    def deriv(T,t):
        T_H1, T_S1, T_H2, T_S2 = T
        dT_H1 = -(Ua + Ub + Uc)*T_H1/Cp_H + Ub*T_H2/Cp_H + Uc*T_S1/Cp_H + P1*u1/Cp_H
        dT_S1 = Uc*T_H1/Cp_S - Uc*T_S1/Cp_S
        dT_H2 = -(Ua + Ub + Uc)*T_H2/Cp_H + Ub*T_H1/Cp_H + Uc*T_S2/Cp_H + P2*u2/Cp_H
        dT_S2 = Uc*T_H2/Cp_S - Uc*T_S2/Cp_S
        return [dT_H1, dT_S1, dT_H2, dT_S2]
    T = odeint(deriv, [0,0,0,0], t)

    plt.figure(figsize=(12,6))
    plt.subplot(2,1,1)
    plt.plot(t, T[:,1] + T_ambient, t, T[:,3] + T_ambient)
    plt.plot(t, T[:,0] + T_ambient)
    plt.plot(t, T[:,2] + T_ambient)
    plt.xlabel('Time / seconds')
    plt.ylabel('Temperature / °C')
    plt.grid()
    plt.subplot(2,1,2)
    plt.plot(t, T[:,1]-T1+T_ambient)
    plt.plot(t, T[:,3]-T2+T_ambient)
    plt.grid
    print(Ua,Ub,Uc,Cp_H,Cp_S)