Temperature Control with an Arduino

This lab is an application of feedback control for a temperature control device. Heat output is adjusted by modulating the voltage to a transistor. A thermistor measures the temperature. Energy from the transistor output is transferred by conduction and convection to the temperature sensor.

This lab teaches principles of system dynamics and control. In particular, this lab reinforces:


Install Python Packages

Step 1: Dynamic Modeling

The three important elements for a control loop are the measurement device (thermistor temperature sensor), an actuator (voltage to the transistor), and capability to perform computerized control (USB interface).

At maximum output the transistor dissipates 3.1 W of power with a voltage of 8.9 V and current of 0.35 A. The transistor is controlled with a 0 to 255 mV signal that linearly adjusts the power output between 0 and 3.1 W. The mass of the transistor and thermistor is about 2 gm with an approximate heat capacity of 460 `J/{kg K}`. The surface area of the heater and sensor is about 5e-4 `m^2`. A convective heat transfer coefficient for quiescent air is approximately 10 `W/{m^2K}`. The heat generated by the transistor transfers primarily by convection and conduction to the temperature sensor. Heat transfer is improved with a thermal coupling (silver epoxy) that connects the two components.

Create a dynamic model of the dynamic response between input power to the transistor and the temperature sensed by the thermistor. Use an energy balance to start the derivation.

$$m\,c_p\frac{dT}{dt} = \sum \dot h_{in} - \sum \dot h_{out} + Q$$

Expand or simplify terms that are needed for this application. The full energy balance includes convection and radiation terms.

$$m\,c_p\frac{dT}{dt} = h\,A\,\left(T_\infty-T\right) + \epsilon\,\sigma\,A\,\left(T_\infty^4-T^4\right) + \alpha\,I\,V$$

where `m` is the mass, `c_p` is the heat capacity, `T` is the temperature, `h` is the heat transfer coefficient, `A` is the area, `T_\infty` is the ambient temperature, `\epsilon` is the emissivity, `\sigma` is the Stefan-Boltzmann constant, `I` is the current, and `V` is the voltage. The parameter `\alpha` is the fraction of electrical energy converted into heat in the transistor. Use this equation to develop a dynamic simulation of the temperature response due to an impulse (off, on, off) in the heater output. Leave the heater on for sufficient time to observe nearly steady state conditions.

Step 2: Generate Data, Fit Physics-based Model

As done in simulation, test the same heater impulse response on the Arduino device and record the temperature. Adjust the uncertain parameters such as `h`, `\alpha`, and `m` from Step 1 to best match the dynamic data from the impulse response. A linearized version of the model can be used to compare the FOPDT model of step 3 to the physics-based model in this step. A linearized model is derived as:

$$m\,c_p\frac{dT}{dt} = f\left(T_\infty,T,V\right) \approx f \left(\bar T_\infty, \bar T, \bar V\right) \ldots \\ + \frac{\partial f}{\partial T_\infty}\bigg|_{\bar T_\infty,\bar T,\bar V} \left(T_\infty-\bar T_\infty\right) \ldots \\ + \frac{\partial f}{\partial T}\bigg|_{\bar T_\infty,\bar T,\bar V} \left(T-\bar T\right) \ldots \\ + \frac{\partial f}{\partial V}\bigg|_{\bar T_\infty,\bar T,\bar V} \left(V-\bar V\right)$$


$$\frac{\partial f}{\partial T_\infty}\bigg|_{\bar T_\infty,\bar T,\bar V} = h\,A\,\bar T_\infty + 4\epsilon\,\sigma\,A\,\bar T_\infty^3 = \beta$$

$$\frac{\partial f}{\partial T}\bigg|_{\bar T_\infty,\bar T,\bar V} = -h\,A\,\bar T - 4\epsilon\,\sigma\,A\,\bar T^3 = \gamma$$

$$\frac{\partial f}{\partial V}\bigg|_{\bar T_\infty,\bar T,\bar V} = \alpha \, \bar I$$

The final linearized equation is further simplified by replacing the partial derivatives with constants `\beta`, `\gamma`, and `\alpha \bar I` that are evaluated at the steady state conditions \bar T_\infty, \bar T, and \bar V.

$$m\,c_p\frac{dT}{dt} = \beta \left(T_\infty-\bar T_\infty\right) + \gamma \left(T-\bar T\right) + \alpha \bar I \left(V-\bar V\right)$$

This is placed into a time-constant form as by dividing by `gamma` and substituting the deviation variables `T\prime_\infty = T_\infty -\bar T_\infty`, `T\prime = T -\bar T`, and `V\prime = V -\bar V`.

$$\frac{m\,c_p}{\gamma}\frac{dT^\prime}{dt} = \frac{\beta}{\gamma} T^\prime_\infty + T^\prime + \frac{\alpha \bar I}{\gamma} V^\prime$$

Step 3: Fit FOPDT Model

Using the data from Step 2, determine the parameters of an FOPDT model that best match the dynamic temperature data including `K_p`, `\tau_p`, and `\theta_p`. A first-order linear system is:

$$\tau_p \frac{dT(t)}{dt} = -T(t) + K_p \, V\left(t-\theta_p\right)$$

where `K_p` is the process gain, `\tau_p` is the process time constant, and `\theta_p` is the process dead-time. A second order (SOPDT) model may also give insight on whether the system shows evidence of a higher order response.

Step 4: PID Control

Obtain PID tuning constants `K_c`, `\tau_I`, and `\tau_D` from IMC correlations. Use the tuning constants for PID control of temperature. Demonstrate step changes in temperature set point and comment on the performance of the controller using the calculated constants. Tune the controller by adjusting the constants to improve performance. Comment on the difference between IMC tuning constants and the improved tuning constants in terms of rise time, overshoot, decay ratio, heater fluctuations, or other relevant performance criteria.

Lab Sample Code

Below is sample code to help with steps of the project. It is a starting point that can be modified to investigate issues such as whether radiative heat transfer is significant, is the temperature response inherently first order or higher order, and what values of uncertain parameters in the physics based model help the predicted temperature agree with the data.

Step 1: Simulate Temperature Response

A first step is to simulate the Arduino temperature response with an energy balance. Just like the physical system, this simulated system responds to changes in the input voltage.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# define Arduino model
def arduino(x,t,mV,Ta):
    ## Parameters
    h = 1.0               # W/m^2-K
    m = 4.0/1000.0        # kg
    Cp = 0.5 * 1000.0     # J/kg-K    
    A = 6.0 / 100.0**2    # Area in m^2
    Amp = 1.0             # Amps to Transitor
    T = x[0]

    # acc = inlet - outlet
    dTdt = (1.0/(m*Cp))*(-h*A*(T-Ta) + Amp*mV/1000.0)
    return dTdt

# Ambient temperature
Ta = (80.0-32.0)*5.0/9.0 + 273.15  # 80 degF in K
# Starting temperature
T0 = Ta  # 80 degF in K

# Time Intervals (sec)
t = np.linspace(0,600,601)

# Manipulated variable
mV = np.ones(len(t))* 0.0 # millivolts input
mV[20:] = 100
# Store results for plotting
T = np.ones(len(t)) * T0

for i in range(len(t)-1):
    ts = [t[i],t[i+1]]
    y = odeint(arduino,T0,ts,args=(mV[i],Ta))
    T0 = y[-1]
    T[i+1] = T0

# Construct results and save data file
# Column 1 = time
# Column 2 = cooling temperature
# Column 3 = reactor temperature
data = np.vstack((t,mV,T)) # vertical stack
data = data.T             # transpose data

# Plot the results

plt.ylabel('Temperature (K)')

plt.xlabel('Time (sec)')

Step 2: Optimize Model Parameters

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import minimize

# Import data file
# Column 1 = time (t)
# Column 2 = input (u)
# Column 3 = output (yp)
data = np.loadtxt('data.txt',delimiter=',')
mV0 = data[0,1]
Tmeas0 = data[0,2]
t = data[:,0].T
mV = data[:,1].T
Tmeas = data[:,2].T

# specify number of steps
ns = len(t)
delta_t = t[1]-t[0]

Cp = 0.5 * 1000.0     # J/kg-K    
A = 6.0 / 100.0**2    # Area in m^2
Ta = Tmeas0           # Ambient temperature (K)

# define Arduino model
def arduino(x,t,mVi,p):
    # uncertain parameters
    h = p[0]
    m = p[1]
    Amp = p[2]

    T = x[0]
    # acc = inlet - outlet
    dTdt = (1.0/(m*Cp))*(-h*A*(T-Ta) + Amp*mVi/1000.0)
    return dTdt

def calc_T(p):
    T = np.ones(len(t)) * Tmeas0
    T0 = T[0]
    for i in range(len(t)-1):
        ts = [t[i],t[i+1]]
        y = odeint(arduino,T0,ts,args=(mV[i],p))
        T0 = y[-1]
        T[i+1] = T0
    return T

# define objective
def objective(p):
    # simulate model
    Tp = calc_T(p)
    # calculate objective
    obj = 0.0
    for i in range(len(Tmeas)):
        obj = obj + ((Tp[i]-Tmeas[i])/Tmeas[i])**2    
    # return result
    return obj

# Parameter initial guess
h = 10.0               # W/m^2-K
m = 4.0/1000.0         # kg
Amp = 1.0              # Amps to transistor
p0 = [h,m,Amp]

# show initial objective
print('Initial SSE Objective: ' + str(objective(p0)))

# optimize parameters
# bounds on variables
bnds = ((1.0, 20.0),(1.0/1000.0,10.0/1000.0),(0.5,2.0))
solution = minimize(objective,p0,method='SLSQP',bounds=bnds)
p = solution.x

# show final objective
print('Final SSE Objective: ' + str(objective(p)))

h = p[0]
m = p[1]
Amp = p[2]
print('h: ' + str(h))
print('m: ' + str(m))
print('Amp: ' + str(Amp))

print('FOPDT Equivalent')
#dx/dt = (-1/taup) * x + (Kp/taup) * u
#dTdt = (1.0/(m*Cp))*(-h*A*(T-Ta) + Amp*mVi/1000.0)
dfdT = -h*A/(m*Cp)
dfdmV = Amp/(1000.0*m*Cp)
taup = -1.0/dfdT
Kp = dfdmV * taup
print('Kp: ' + str(Kp))
print('taup: ' + str(taup))

# calculate model with updated parameters
T1 = calc_T(p0)
T2 = calc_T(p)

# Plot the results

plt.plot(t,T1,'b:',linewidth=3,label='Initial Guess')
plt.plot(t,T2,'k--',linewidth=3,label='Final Prediction')
plt.ylabel('Temperature (K)')
plt.xlabel('Time (sec)')

Step 3: Fit FOPDT and SOPDT to Data

Determine the parameters of an FOPDT model that best match the dynamic temperature data including `K_p`, `\tau_p`, and `\theta_p`. A second order (SOPDT) can also be fit to investigate whether a higher order model is more accurate.

Step 4: Implement PID Controller

Simulate a PID Controller with one of the models determined in steps 2 or 3. It is suggested to tune the controller in simulation before implementing with an Arduino. Tuning on a device that takes 10-20 minutes per test is much slower than running a PID controller in simulation. Once optimized PID tuning values are obtained, demonstrate the performance with the physical device.

Solution Files

Additional approaches for modeling, estimation, and control are provided at the link below. These include:

  • Dynamic modeling
    • First principles
    • Empirical / machine learning
  • Dynamic Estimation
    • Parameter regression
    • Kalman filter
    • Moving Horizon Estimation
  • Process control
    • ON/OFF control
    • PID control
    • Stability analysis
  • Advanced control
    • Linear and Nonlinear
    • Model Predictive Control
    • SISO - Single Input Single Output
    • MIMO - Multiple Input Multiple Output