The objective of this activity is to fit the physics-based predictions to the data as well as fit a first-order plus dead-time (FOPDT) model to the data. In both cases, select parameters are adjusted to minimize an objective function such as a sum of squared errors between the model predicted values and the measured values. Test the temperature response of the Arduino device by introducing a step in the heater.
Fit FOPDT Model with Optimization
Determine the parameters of an FOPDT model that best match the dynamic temperature data including KpKp, τpτp, and θpθp. A second order (SOPDT) can also be fit to investigate whether a higher order model is more accurate.
import numpy as np import matplotlib.pyplotas plt from scipy.integrateimport odeint from scipy.optimizeimport minimize from scipy.interpolateimport interp1d
# Import CSV data file # Column 1 = time (t) # Column 2 = input (u) # Column 3 = output (yp)
data = np.loadtxt('data.txt',delimiter=',',skiprows=1)
u0 = data[0,1]
yp0 = data[0,3]
t = data[:,0].T - data[0,0]
u = data[:,1].T
yp = data[:,3].T
# specify number of steps
ns =len(t)
delta_t = t[1]-t[0] # create linear interpolation of the u data versus time
uf = interp1d(t,u)
# define first-order plus dead-time approximation def fopdt(y,t,uf,Km,taum,thetam): # arguments # y = output # t = time # uf = input linear function (for time shift) # Km = model gain # taum = model time constant # thetam = model time constant # time-shift u try: if(t-thetam)<=0:
um = uf(0.0) else:
um = uf(t-thetam) except: #print('Error with time extrapolation: ' + str(t))
um = u0 # calculate derivative
dydt =(-(y-yp0) + Km * (um-u0))/taum return dydt
# simulate FOPDT model with x=[Km,taum,thetam] def sim_model(x): # input arguments
Km = x[0]
taum = x[1]
thetam = x[2] # storage for model values
ym = np.zeros(ns)# model # initial condition
ym[0]= yp0 # loop through time steps for i inrange(0,ns-1):
ts =[t[i],t[i+1]]
y1 = odeint(fopdt,ym[i],ts,args=(uf,Km,taum,thetam))
ym[i+1]= y1[-1] return ym
# define objective def objective(x): # simulate model
ym = sim_model(x) # calculate objective
obj =0.0 for i inrange(len(ym)):
obj = obj + (ym[i]-yp[i])**2 # return result return obj
# show initial objective print('Initial SSE Objective: ' + str(objective(x0)))
# optimize Km, taum, thetam
solution = minimize(objective,x0)
# Another way to solve: with bounds on variables #bnds = ((0.4, 0.6), (1.0, 10.0), (0.0, 30.0)) #solution = minimize(objective,x0,bounds=bnds,method='SLSQP')
x = solution.x
# show final objective print('Final SSE Objective: ' + str(objective(x)))
where mm is the mass, cpcp is the heat capacity, TT is the temperature, UU is the heat transfer coefficient, AA is the area, T∞T∞ is the ambient temperature, ε=0.9ε=0.9 is the emissivity, σ=σ= 5.67x10-8Wm2K4Wm2K4 is the Stefan-Boltzmann constant, and QQ is the percentage heater output. The parameter αα is a factor that relates heater output (0-100%) to power dissipated by the transistor in Watts.
Adjust the uncertain parameters such as UU and αα from the modeling exercise to best match the dynamic data from the impulse response.
import numpy as np import matplotlib.pyplotas plt from scipy.integrateimport odeint from scipy.optimizeimport minimize
# Import data file # Column 1 = time (t) # Column 2 = input (u) # Column 3 = output (yp)
data = np.loadtxt('data.txt',delimiter=',',skiprows=1) # initial conditions
Q0 = data[0,1]
Tmeas0 = data[0,3] # extract data columns
t = data[:,0].T
Q = data[:,1].T
Tmeas = data[:,3].T
# specify number of steps
ns =len(t)
delta_t = t[1]-t[0]
Cp =0.5 * 1000.0# J/kg-K
A =12.0 / 100.0**2# Area in m^2
Ta = Tmeas0 # Ambient temperature (K)
# define energy balance model def heat(x,t,Q,p): # Adjustable Parameters
U = p[0]# starting at 10.0 W/m^2-K
alpha = p[1]# starting as 0.01 W / % heater
# Known Parameters
m =4.0/1000.0# kg
Cp =0.5 * 1000.0# J/kg-K
A =12.0 / 100.0**2# Area in m^2
eps =0.9# Emissivity
sigma =5.67e-8# Stefan-Boltzman
def calc_T(p):
T = np.ones(len(t)) * Tmeas0
T0 = T[0] for i inrange(len(t)-1):
ts =[t[i],t[i+1]]
y = odeint(heat,T0,ts,args=(Q[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 inrange(len(Tmeas)):
obj = obj + ((Tp[i]-Tmeas[i])/Tmeas[i])**2 # return result return obj
# Parameter initial guess
U =10.0# Heat transfer coefficient (W/m^2-K)
alpha =0.01# Heat gain (W/%)
p0 =[U,alpha]
# show initial objective print('Initial SSE Objective: ' + str(objective(p0)))
# optimize parameters # bounds on variables
bnds =((2.0,20.0),(0.005,0.02))
solution = minimize(objective,p0,method='SLSQP',bounds=bnds)
p = solution.x
# show final objective print('Final SSE Objective: ' + str(objective(p)))
# Known Parameters
m =4.0/1000.0# kg
Cp =0.5 * 1000.0# J/kg-K
A =12.0 / 100.0**2# Area in m^2
eps =0.9# Emissivity
sigma =5.67e-8# Stefan-Boltzman
T0 = Tmeas0
The final linearized equation is further simplified by replacing the partial derivatives with constants ββ, γγ, and αα that are evaluated at the steady state conditions ˉT∞¯¯¯T∞, ˉT¯¯¯T, and ˉQ¯¯¯Q.
This is placed into a time-constant form as by dividing by γγ and substituting the deviation variables T′∞=T∞-ˉT∞, T′=T-ˉT, and Q′=Q-ˉQ.
mcpγdT′dt=βγT′∞+T′+αγQ′
Further simplification is possible with β=-γ.
mcpγdT′dt=−T′∞+T′+αγQ′
Multiplying both sides by -1 puts the equation into a time constant form with new constants τP and Kp.
τPdT′dt=−T′+T′∞+KpQ′
KP=−αγτP=−mcpγ
This derivation does not include a time delay but the other parameters can be correlated to an empirical FOPDT fit of the data. Another way to fit the uncertain parameters such as overall heat transfer coefficient and α is to use optimization methods. There is a script that demonstrates the use of optimization to minimize an objective function. The objective function is a minimization of the difference between the predicted and measured values.