Main

## Parameter Estimation with Arduino Data

The objective of this activity is to fit the physics-based predictions to the data for a two heater model of the temperature control lab. 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. Use an optimizer to adjust the parameters and achieve alignment between the model and the measured values.

Step Test to Obtain Data

Start the exercise by collecting heater step test data from the device.

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

# define energy balance model
def heat(x,t,Q1,Q2):
# 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 = 10.0 / 100.0**2 # Area in m^2
As = 2.0 / 100.0**2 # Area in m^2
alpha1 = 0.0100     # W / % heater 1
alpha2 = 0.0075     # W / % heater 2
eps = 0.9          # Emissivity
sigma = 5.67e-8    # Stefan-Boltzman

# Temperature States
T1 = x[0]
T2 = x[1]

# Heat Transfer Exchange Between 1 and 2
conv12 = U*As*(T2-T1)
rad12  = eps*sigma*As * (T2**4 - T1**4)

# Nonlinear Energy Balances
dT1dt = (1.0/(m*Cp))*(U*A*(Ta-T1) \
+ eps * sigma * A * (Ta**4 - T1**4) \
+ conv12 + rad12 \
+ alpha1*Q1)
dT2dt = (1.0/(m*Cp))*(U*A*(Ta-T2) \
+ eps * sigma * A * (Ta**4 - T2**4) \
- conv12 - rad12 \
+ alpha2*Q2)

return [dT1dt,dT2dt]

# save txt file
def save_txt(t,u1,u2,y1,y2,sp1,sp2):
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)'

# Connect to Arduino
a = tclab.TCLab()

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

# Run time in minutes
run_time = 10.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
Tp1 = np.ones(loops) * a.T1
Tp2 = np.ones(loops) * a.T2
error_eb = np.zeros(loops)

# impulse tests (0 - 100%)
Q1 = np.ones(loops) * 0.0
Q2 = np.ones(loops) * 0.0

# steps
Q1[10:] = 100.0
Q2[100:] = 50.0
Q1[200:] = 5.0
Q2[300:] = 80.0
Q1[400:] = 70.0
Q2[500:] = 10.0

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
Tinit = [Tp1[i-1]+273.15,Tp2[i-1]+273.15]
Tnext = odeint(heat,Tinit, \
[0,dt],args=(Q1[i-1],Q2[i-1]))
Tp1[i] = Tnext[1,0]-273.15
Tp2[i] = Tnext[1,1]-273.15
error_eb[i] = error_eb[i-1] \
+ (abs(Tp1[i]-T1[i]) \
+  abs(Tp2[i]-T2[i]))*dt

# 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],Tp1[0:i],'k-',label=r'$T_1$ energy balance')
plt.plot(tm[0:i],T2[0:i],'bx',label=r'$T_2$ measured')
plt.plot(tm[0:i],Tp2[0:i],'k--',label=r'$T_2$ energy balance')
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 Error')
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 and plot at end
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.png')

# 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.png')

# 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.png')
raise

Fit Dual Heater Model with Optimization

Starting with the dual heater energy balances, adjust parameters such as U the heat transfer coefficient, \alpha_1 or \alpha_2 the heater power factors, or other uncertain parameters that will improve the fit to the data. Other factors such as the ambient temperature (T_\infty), can be updated by noting the temperature when the device is cool and the heaters are initially off.

Parameter Estimation with Python Scipy

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

# generate data file from TCLab or get sample data file from:
#  https://apmonitor.com/pdc/index.php/Main/ArduinoEstimation2
# Import data file
# Column 1 = time (t)
# Column 2 = input (u)
# Column 3 = output (yp)
# extract data columns
t = data[:,0].T
Q1 = data[:,1].T
Q2 = data[:,2].T
T1meas = data[:,3].T
T2meas = data[:,4].T

# number of time points
ns = len(t)

# define energy balance model
def heat(x,t,Q1,Q2,p):
# Optimized parameters
U,alpha1,alpha2 = p

# Parameters
Ta = 23 + 273.15   # K
m = 4.0/1000.0     # kg
Cp = 0.5 * 1000.0  # J/kg-K
A = 10.0 / 100.0**2 # Area in m^2
As = 2.0 / 100.0**2 # Area in m^2
eps = 0.9          # Emissivity
sigma = 5.67e-8    # Stefan-Boltzman

# Temperature States
T1 = x[0] + 273.15
T2 = x[1] + 273.15

# Heat Transfer Exchange Between 1 and 2
conv12 = U*As*(T2-T1)
rad12  = eps*sigma*As * (T2**4 - T1**4)

# Nonlinear Energy Balances
dT1dt = (1.0/(m*Cp))*(U*A*(Ta-T1) \
+ eps * sigma * A * (Ta**4 - T1**4) \
+ conv12 + rad12 \
+ alpha1*Q1)
dT2dt = (1.0/(m*Cp))*(U*A*(Ta-T2) \
+ eps * sigma * A * (Ta**4 - T2**4) \
- conv12 - rad12 \
+ alpha2*Q2)

return [dT1dt,dT2dt]

def simulate(p):
T = np.zeros((len(t),2))
T[0,0] = T1meas[0]
T[0,1] = T2meas[0]
T0 = T[0]
for i in range(len(t)-1):
ts = [t[i],t[i+1]]
y = odeint(heat,T0,ts,args=(Q1[i],Q2[i],p))
T0 = y[-1]
T[i+1] = T0
return T

# define objective
def objective(p):
# simulate model
Tp = simulate(p)
# calculate objective
obj = 0.0
for i in range(len(t)):
obj += ((Tp[i,0]-T1meas[i])/T1meas[i])**2 \
+((Tp[i,1]-T2meas[i])/T2meas[i])**2
# return result
return obj

# Parameter initial guess
U = 10.0           # Heat transfer coefficient (W/m^2-K)
alpha1 = 0.0100    # Heat gain 1 (W/%)
alpha2 = 0.0075    # Heat gain 2 (W/%)
p0 = [U,alpha1,alpha2]

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

# optimize parameters
# bounds on variables
bnds = ((2.0, 20.0),(0.005,0.02),(0.002,0.015))
solution = minimize(objective,p0,method='SLSQP',bounds=bnds)
p = solution.x

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

# optimized parameter values
U = p[0]
alpha1 = p[1]
alpha2 = p[2]
print('U: ' + str(U))
print('alpha1: ' + str(alpha1))
print('alpha2: ' + str(alpha2))

# calculate model with updated parameters
Ti  = simulate(p0)
Tp  = simulate(p)

# Plot results
plt.figure(1)

plt.subplot(3,1,1)
plt.plot(t/60.0,Ti[:,0],'y:',label=r'$T_1$ initial')
plt.plot(t/60.0,T1meas,'b-',label=r'$T_1$ measured')
plt.plot(t/60.0,Tp[:,0],'r--',label=r'$T_1$ optimized')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')

plt.subplot(3,1,2)
plt.plot(t/60.0,Ti[:,1],'y:',label=r'$T_2$ initial')
plt.plot(t/60.0,T2meas,'b-',label=r'$T_2$ measured')
plt.plot(t/60.0,Tp[:,1],'r--',label=r'$T_2$ optimized')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')

plt.subplot(3,1,3)
plt.plot(t/60.0,Q1,'g-',label=r'$Q_1$')
plt.plot(t/60.0,Q2,'k--',label=r'$Q_2$')
plt.ylabel('Heater Output')
plt.legend(loc='best')

plt.xlabel('Time (min)')
plt.show()

Parameter Estimation with Python GEKKO

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

# Import or generate data
filename = 'tclab_dyn_data2.csv'
try:
except:

# Create GEKKO Model
m = GEKKO()
m.time = data['Time'].values

# Parameters to Estimate
U = m.FV(value=10,lb=1,ub=20)
alpha1 = m.FV(value=0.01,lb=0.003,ub=0.03)  # W / % heater
alpha2 = m.FV(value=0.005,lb=0.002,ub=0.02) # W / % heater

# STATUS=1 allows solver to adjust parameter
U.STATUS = 1
alpha1.STATUS = 1
alpha2.STATUS = 1

# Measured inputs
Q1 = m.MV(value=data['H1'].values)
Q2 = m.MV(value=data['H2'].values)

# State variables
TC1 = m.CV(value=data['T1'].values)
TC1.FSTATUS = 1    # minimize fstatus * (meas-pred)^2
TC2 = m.CV(value=data['T2'].values)
TC2.FSTATUS = 1    # minimize fstatus * (meas-pred)^2

Ta = m.Param(value=19.0+273.15)     # K
mass = m.Param(value=4.0/1000.0)    # kg
Cp = m.Param(value=0.5*1000.0)      # J/kg-K
A = m.Param(value=10.0/100.0**2)    # Area not between heaters in m^2
As = m.Param(value=2.0/100.0**2)    # Area between heaters in m^2
eps = m.Param(value=0.9)            # Emissivity
sigma = m.Const(5.67e-8)            # Stefan-Boltzmann

# Heater temperatures in Kelvin
T1 = m.Intermediate(TC1+273.15)
T2 = m.Intermediate(TC2+273.15)

# Heat transfer between two heaters
Q_C12 = m.Intermediate(U*As*(T2-T1)) # Convective
Q_R12 = m.Intermediate(eps*sigma*As*(T2**4-T1**4)) # Radiative

# Semi-fundamental correlations (energy balances)
m.Equation(TC1.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T1) \
+ eps * sigma * A * (Ta**4 - T1**4) \
+ Q_C12 + Q_R12 \
+ alpha1*Q1))

m.Equation(TC2.dt() == (1.0/(mass*Cp))*(U*A*(Ta-T2) \
+ eps * sigma * A * (Ta**4 - T2**4) \
- Q_C12 - Q_R12 \
+ alpha2*Q2))

# Options
m.options.IMODE   = 5 # MHE
m.options.EV_TYPE = 2 # Objective type
m.options.NODES   = 2 # Collocation nodes
m.options.SOLVER  = 3 # IPOPT

# Solve
m.solve(disp=True)

# Parameter values
print('U     : ' + str(U.value[0]))
print('alpha1: ' + str(alpha1.value[0]))
print('alpha2: ' + str(alpha2.value[0]))

# Create plot
plt.figure()
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(data['Time'],data['T1'],'ro',label=r'$T_1$ measured')
plt.plot(m.time,TC1.value,color='purple',linestyle='--',\
linewidth=3,label=r'$T_1$ predicted')
plt.plot(data['Time'],data['T2'],'bx',label=r'$T_2$ measured')
plt.plot(m.time,TC2.value,color='orange',linestyle='--',\
linewidth=3,label=r'$T_2$ predicted')
plt.ylabel('Temperature (degC)')
plt.legend(loc=2)
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(data['Time'],data['H1'],'r-',\
linewidth=3,label=r'$Q_1$')
plt.plot(data['Time'],data['H2'],'b:',\
linewidth=3,label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.show()

Parameter Estimation with Confidence Intervals

See Regression Statistics Introduction for information on obtaining parameter confidence intervals and data prediction bands.

import numpy as np
from scipy.optimize import curve_fit
import uncertainties as unc
import matplotlib.pyplot as plt
import uncertainties.unumpy as unp
from scipy.integrate import odeint
import pandas as pd
from scipy import stats

# calculate lower and upper prediction bands
def predband(x, xd, yd, f_vars, conf=0.95):
"""
Code adapted from Rodrigo Nemmen's post:
https://astropython.blogspot.com.ar/2011/12/calculating-prediction-band-
of-linear.html

Calculates the prediction band of the regression model at the
desired confidence level.

Clarification of the difference between confidence and prediction bands:

"The prediction bands are further from the best-fit line than the
confidence bands, a lot further if you have many data points. The 95%
prediction band is the area in which you expect 95% of all data points
to fall. In contrast, the 95% confidence band is the area that has a
95% chance of containing the true regression line."
reg_graphing_tips_linear_regressio.htm)

Arguments:
- x: array with x values to calculate the confidence band.
- xd, yd: data arrays.
- a, b, c: linear fit parameters.
- conf: desired confidence level, by default 0.95 (2 sigma)

References:
1. https://www.JerryDallal.com/LHSP/slr.htm, Introduction to Simple Linear
Regression, Gerard E. Dallal, Ph.D.
"""

alpha = 1. - conf    # Significance
N = xd.size          # data sample size
var_n = len(f_vars)  # Number of variables used by the fitted function.
# Quantile of Student's t distribution for p=(1 - alpha/2)
q = stats.t.ppf(1. - alpha / 2., N - var_n)
# Std. deviation of an individual measurement (Bevington, eq. 6.15)
se = np.sqrt(1. / (N - var_n) * np.sum((yd - simulate(xd, *f_vars)) ** 2))
# Auxiliary definitions
sx = (x - xd.mean()) ** 2
sxd = np.sum((xd - xd.mean()) ** 2)
# Predicted values (best-fit model)
yp = simulate(x, *f_vars)
# Prediction band
dy = q * se * np.sqrt(1. + (1. / N) + (sx / sxd))
# Upper & lower prediction bands.
lpb, upb = yp - dy, yp + dy
return lpb, upb

# generate data file from TCLab or get sample data file from:
#  https://apmonitor.com/pdc/index.php/Main/ArduinoEstimation2
# Import data file
# Column 1 = time (t)
# Column 2 = input (u)
# Column 3 = output (yp)
# extract data columns
t = data[:,0].T
Q1 = data[:,1].T
Q2 = data[:,2].T
T1meas = data[:,3].T
T2meas = data[:,4].T
ind = np.linspace(0,np.size(t),np.size(t))

# number of time points
ns = len(t)

# define energy balance model
def heat(x,t,Q1,Q2,p):
# Optimized parameters
U,alpha1,alpha2 = p

# Parameters
Ta = 23 + 273.15   # K
m = 4.0/1000.0     # kg
Cp = 0.5 * 1000.0  # J/kg-K
A = 10.0 / 100.0**2 # Area in m^2
As = 2.0 / 100.0**2 # Area in m^2
eps = 0.9          # Emissivity
sigma = 5.67e-8    # Stefan-Boltzman

# Temperature States
T1 = x[0] + 273.15
T2 = x[1] + 273.15

# Heat Transfer Exchange Between 1 and 2
conv12 = U*As*(T2-T1)
rad12  = eps*sigma*As * (T2**4 - T1**4)

# Nonlinear Energy Balances
dT1dt = (1.0/(m*Cp))*(U*A*(Ta-T1) \
+ eps * sigma * A * (Ta**4 - T1**4) \
+ conv12 + rad12 \
+ alpha1*Q1)
dT2dt = (1.0/(m*Cp))*(U*A*(Ta-T2) \
+ eps * sigma * A * (Ta**4 - T2**4) \
- conv12 - rad12 \
+ alpha2*Q2)

return [dT1dt,dT2dt]

def simulate(tm,U,alpha1,alpha2):
T = np.zeros((len(t),2))
T[0,0] = T1meas[0]
T[0,1] = T2meas[0]
T0 = T[0]
p = (U,alpha1,alpha2)
for i in range(len(t)-1):
ts = [t[i],t[i+1]]
y = odeint(heat,T0,ts,args=(Q1[i],Q2[i],p))
T0 = y[-1]
T[i+1] = T0
z = np.empty((len(t)*2))
z[0:len(t)] = T[:,0]
z[len(t):] = T[:,1]
return z

def simulate2(p):
T = np.zeros((len(t),2))
T[0,0] = T1meas[0]
T[0,1] = T2meas[0]
T0 = T[0]
for i in range(len(t)-1):
ts = [t[i],t[i+1]]
y = odeint(heat,T0,ts,args=(Q1[i],Q2[i],p))
T0 = y[-1]
T[i+1] = T0
return T

# Parameter initial guess
U = 10.0           # Heat transfer coefficient (W/m^2-K)
alpha1 = 0.0100    # Heat gain 1 (W/%)
alpha2 = 0.0075    # Heat gain 2 (W/%)
pinit = [U,alpha1,alpha2]

x = []
y = np.empty((len(t)*2))
y[0:len(t)] = T1meas
y[len(t):] = T2meas

popt, pcov = curve_fit(simulate, x, y)

Uu, alpha1u, alpha2u = unc.correlated_values(popt, pcov)

# create prediction band
lpb, upb = predband(y, y, y, popt, conf=0.95)
lpb1 = np.empty((len(t)))
lpb2 = np.empty((len(t)))
upb1 = np.empty((len(t)))
upb2 = np.empty((len(t)))
lpb1[0:len(t)] = lpb[0:len(t)]
lpb2[0:len(t)] = lpb[len(t):]
upb1[0:len(t)] = upb[0:len(t)]
upb2[0:len(t)] = upb[len(t):]

# optimized parameter values with uncertainties
print('Optimal Parameters with Uncertanty Range')
print('U: ' + str(Uu))
print('alpha1: ' + str(alpha1u))
print('alpha2: ' + str(alpha2u))

# calculate model with updated parameters
Ti  = simulate2(pinit)
Tp  = simulate2(popt)

# Plot results
plt.figure(1)

plt.subplot(3,1,1)
plt.plot(t/60.0,Ti[:,0],'y:',label=r'$T_1$ initial')
plt.plot(t/60.0,T1meas,'b-',label=r'$T_1$ measured')
plt.plot(t/60.0,Tp[:,0],'r--',label=r'$T_1$ optimized')
plt.plot(t/60.0,lpb1,'k:',label=r'$T_1$ prediction band')
plt.plot(t/60.0,upb1,'k:')

plt.ylabel('Temperature (degC)')
plt.legend(loc='best')

plt.subplot(3,1,2)
plt.plot(t/60.0,Ti[:,1],'y:',label=r'$T_2$ initial')
plt.plot(t/60.0,T2meas,'b-',label=r'$T_2$ measured')
plt.plot(t/60.0,Tp[:,1],'r--',label=r'$T_2$ optimized')
plt.plot(t/60.0,lpb2,'k:',label=r'$T_2$ prediction band')
plt.plot(t/60.0,upb2,'k:')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')

plt.subplot(3,1,3)
plt.plot(t/60.0,Q1,'g-',label=r'$Q_1$')
plt.plot(t/60.0,Q2,'k--',label=r'$Q_2$')
plt.ylabel('Heater Output')
plt.legend(loc='best')

plt.xlabel('Time (min)')
plt.show()

Parameter Estimation with MATLAB fmincon

close all; clear all; clc

global t T1meas T2meas Q1 Q2

% generate data file from TCLab or get sample data file from:
%  https://apmonitor.com/pdc/index.php/Main/ArduinoEstimation2
% Import data file
% Column 1 = time (t)
% Column 2 = input (u)
% Column 3 = output (yp)
filename = 'data.txt';
delimiterIn = ',';
% extract data columns
t = z.data(:,1);
Q1 = z.data(:,2);
Q2 = z.data(:,3);
T1meas = z.data(:,4);
T2meas = z.data(:,5);

% number of time points
ns = length(t);

% Parameter initial guess
U = 10.0;           % Heat transfer coefficient (W/m^2-K)
alpha1 = 0.0100;    % Heat gain 1 (W/%)
alpha2 = 0.0075;    % Heat gain 2 (W/%)
p0 = [U,alpha1,alpha2];

% show initial objective
disp(['Initial SSE Objective: ' num2str(objective(p0))])

% optimize parameters

% no linear constraints
A = [];
b = [];
Aeq = [];
beq = [];
nlcon = [];

% optimize with fmincon
lb = [2,0.005,0.002]; % lower bound
ub = [20,0.02,0.015]; % upper bound
% options = optimoptions(@fmincon,'Algorithm','interior-point');
p = fmincon(@objective,p0,A,b,Aeq,beq,lb,ub,nlcon); %,options);

% show final objective
disp(['Final SSE Objective: ' num2str(objective(p))])

% optimized parameter values
U = p(1);
alpha1 = p(2);
alpha2 = p(3);
disp(['U: ' num2str(U)])
disp(['alpha1: ' num2str(alpha1)])
disp(['alpha2: ' num2str(alpha2)])

% calculate model with updated parameters
Ti  = simulate(p0);
Tp  = simulate(p);

% Plot results
figure(1)

subplot(3,1,1)
plot(t/60.0,Ti(:,1),'y:','LineWidth',2)
hold on
plot(t/60.0,T1meas,'b-','LineWidth',2)
plot(t/60.0,Tp(:,1),'r--','LineWidth',2)
ylabel('Temperature (degC)')
legend('T_1 initial','T_1 measured','T_1 optimized')

subplot(3,1,2)
plot(t/60.0,Ti(:,2),'y:','LineWidth',2)
hold on
plot(t/60.0,T2meas,'b-','LineWidth',2)
plot(t/60.0,Tp(:,2),'r--','LineWidth',2)
ylabel('Temperature (degC)')
legend('T_2 initial','T_2 measured','T_2 optimized')

subplot(3,1,3)
plot(t/60.0,Q1,'g-','LineWidth',2)
hold on
plot(t/60.0,Q2,'k--','LineWidth',2)
ylabel('Heater Output')
legend('Q_1','Q_2')

xlabel('Time (min)')

Save as heat.m

% define energy balance model
function dTdt = heat(t,x,Q1,Q2,p)
%U = 10.0;           % W/m^2-K
%alpha1 = 0.0100;    % W / % heater 1
%alpha2 = 0.0075;    % W / % heater 2
U = p(1);
alpha1 = p(2);
alpha2 = p(3);

% Parameters
Ta = 23 + 273.15;   % K
m = 4.0/1000.0;     % kg
Cp = 0.5 * 1000.0;  % J/kg-K
A = 10.0 / 100.0^2; % Area in m^2
As = 2.0 / 100.0^2; % Area in m^2
eps = 0.9;          % Emissivity
sigma = 5.67e-8;    % Stefan-Boltzman

% Temperature States
T1 = x(1)+273.15;
T2 = x(2)+273.15;

% Heat Transfer Exchange Between 1 and 2
conv12 = U*As*(T2-T1);
rad12  = eps*sigma*As * (T2^4 - T1^4);

% Nonlinear Energy Balances
dT1dt = (1.0/(m*Cp))*(U*A*(Ta-T1) ...
+ eps * sigma * A * (Ta^4 - T1^4) ...
+ conv12 + rad12 ...
+ alpha1*Q1);
dT2dt = (1.0/(m*Cp))*(U*A*(Ta-T2) ...
+ eps * sigma * A * (Ta^4 - T2^4) ...
- conv12 - rad12 ...
+ alpha2*Q2);
dTdt = [dT1dt,dT2dt]';
end

Save as simulate.m

function T = simulate(p)

global t T1meas T2meas Q1 Q2

T = zeros(length(t),2);
T(1,1) = T1meas(1);
T(1,2) = T2meas(1);
T0 = [T(1,1),T(1,2)];
for i = 1:length(t)-1
ts = [t(i),t(i+1)];
sol = ode15s(@(t,x)heat(t,x,Q1(i),Q2(i),p),ts,T0);
T0 = [sol.y(1,end),sol.y(2,end)];
T(i+1,1) = T0(1);
T(i+1,2) = T0(2);
end

end

Save as objective.m

% define objective

function obj = objective(p)
global T1meas T2meas
% simulate model
Tp = simulate(p);
% calculate objective
obj =  sum(((Tp(:,1)-T1meas)./T1meas).^2) ...
+ sum(((Tp(:,2)-T2meas)./T2meas).^2);
end