import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from gekko import GEKKO
import time
# -------------------------------------
# import or generate data
# -------------------------------------
filename = 'tclab_ss_data2.csv'
try:
try:
data = pd.read_csv(filename)
except:
url = 'https://apmonitor.com/do/uploads/Main/tclab_ss_data2.txt'
data = pd.read_csv(url)
except:
# generate training data if data file not available
import tclab
# Connect to Arduino
a = tclab.TCLab()
fid = open(filename,'w')
fid.write('Heater 1,Heater 2,Temperature 1,Temperature 2\n')
fid.close()
# data collection takes 6 hours = 120 pts * 3 minutes each
npts = 120
for i in range(npts):
# set random heater values
Q1 = np.random.rand()*100
Q2 = np.random.rand()*100
a.Q1(Q1)
a.Q2(Q2)
print('Heater 1: ' + str(Q1) + ' %')
print('Heater 2: ' + str(Q2) + ' %')
# wait 3 minutes
time.sleep(3*60)
# record temperature and heater value
print('Temperature 1: ' + str(a.T1) + ' degC')
print('Temperature 2: ' + str(a.T2) + ' degC')
fid = open(filename,'a')
fid.write(str(Q1)+','+str(Q2)+','+str(a.T1)+','+str(a.T2)+'\n')
fid.close()
# close connection to Arduino
a.close()
# read data file
data = pd.read_csv(filename)
# -------------------------------------
# scale data
# -------------------------------------
s = MinMaxScaler(feature_range=(0,1))
sc_train = s.fit_transform(data)
# partition into inputs and outputs
xs = sc_train[:,0:2] # 2 heaters
ys = sc_train[:,2:4] # 2 temperatures
# -------------------------------------
# build neural network
# -------------------------------------
nin = 2 # inputs
n1 = 2 # hidden layer 1 (linear)
n2 = 2 # hidden layer 2 (nonlinear)
n3 = 2 # hidden layer 3 (linear)
nout = 2 # outputs
# Initialize gekko models
train = GEKKO()
dyn = GEKKO()
model = [train,dyn]
for m in model:
# use APOPT solver
m.options.SOLVER = 1
# input(s)
m.inpt = [m.Param() for i in range(nin)]
# layer 1 (linear)
m.w1 = m.Array(m.FV, (nout,nin,n1))
m.l1 = [[m.Intermediate(sum([m.w1[k,j,i]*m.inpt[j] \
for j in range(nin)])) for i in range(n1)] \
for k in range(nout)]
# layer 2 (tanh)
m.w2 = m.Array(m.FV, (nout,n1,n2))
m.l2 = [[m.Intermediate(sum([m.tanh(m.w2[k,j,i]*m.l1[k][j]) \
for j in range(n1)])) for i in range(n2)] \
for k in range(nout)]
# layer 3 (linear)
m.w3 = m.Array(m.FV, (nout,n2,n3))
m.l3 = [[m.Intermediate(sum([m.w3[k,j,i]*m.l2[k][j] \
for j in range(n2)])) for i in range(n3)] \
for k in range(nout)]
# outputs
m.outpt = [m.CV() for i in range(nout)]
m.Equations([m.outpt[k]==sum([m.l3[k][i] for i in range(n3)]) \
for k in range(nout)])
# flatten matrices
m.w1 = m.w1.flatten()
m.w2 = m.w2.flatten()
m.w3 = m.w3.flatten()
# -------------------------------------
# fit parameter weights
# -------------------------------------
m = train
for i in range(nin):
m.inpt[i].value=xs[:,i]
for i in range(nout):
m.outpt[i].value = ys[:,i]
m.outpt[i].FSTATUS = 1
for i in range(len(m.w1)):
m.w1[i].FSTATUS=1
m.w1[i].STATUS=1
m.w1[i].MEAS=1.0
for i in range(len(m.w2)):
m.w2[i].STATUS=1
m.w2[i].FSTATUS=1
m.w2[i].MEAS=0.5
for i in range(len(m.w3)):
m.w3[i].FSTATUS=1
m.w3[i].STATUS=1
m.w3[i].MEAS=1.0
m.options.IMODE = 2
m.options.EV_TYPE = 2
# solve for weights to minimize loss (objective)
m.solve(disp=True)
# -------------------------------------
# generate dynamic predictions
# -------------------------------------
m = dyn
tf = 600
m.time = np.linspace(0,tf,tf+1)
# load neural network parameters
for i in range(len(m.w1)):
m.w1[i].MEAS=train.w1[i].NEWVAL
m.w1[i].FSTATUS = 1
for i in range(len(m.w2)):
m.w2[i].MEAS=train.w2[i].NEWVAL
m.w2[i].FSTATUS = 1
for i in range(len(m.w3)):
m.w3[i].MEAS=train.w3[i].NEWVAL
m.w3[i].FSTATUS = 1
# step tests
Q1d = np.zeros(tf+1)
Q1d[10:200] = 80
Q1d[200:280] = 20
Q1d[280:400] = 70
Q1d[400:] = 50
Q1 = m.Param()
Q1.value = Q1d
Q2d = np.zeros(tf+1)
Q2d[120:320] = 100
Q2d[320:520] = 10
Q2d[520:] = 80
Q2 = m.Param()
Q2.value = Q2d
# scaled inputs
m.inpt[0].value = Q1d * s.scale_[0] + s.min_[0]
m.inpt[1].value = Q2d * s.scale_[1] + s.min_[1]
# define Temperature output
Q0 = 0 # initial heater
T0 = 19 # ambient temperature
# scaled steady state ouput
T1_ss = m.Var(value=T0)
T2_ss = m.Var(value=T0)
m.Equation(T1_ss == (m.outpt[0]-s.min_[2])/s.scale_[2])
m.Equation(T2_ss == (m.outpt[1]-s.min_[3])/s.scale_[3])
# dynamic prediction
T1 = m.Var(value=T0)
T2 = m.Var(value=T0)
# time constant
tau = m.Param(value=120) # determine in a later exercise
# additional model equation for dynamics
m.Equation(tau*T1.dt()==-(T1-T0)+(T1_ss-T0))
m.Equation(tau*T2.dt()==-(T2-T0)+(T2_ss-T0))
# solve dynamic simulation
m.options.IMODE=4
m.solve()
# generate step test data on Arduino
# -------------------------------------
# import or generate data
# -------------------------------------
filename = 'tclab_dyn_data2.csv'
try:
try:
data = pd.read_csv(filename)
except:
url = 'https://apmonitor.com/do/uploads/Main/tclab_dyn_data2.txt'
data = pd.read_csv(url)
except:
# generate training data if data file not available
import tclab
# Connect to Arduino
a = tclab.TCLab()
fid = open(filename,'w')
fid.write('Time,H1,H2,T1,T2\n')
fid.close()
# check for cool down
i = 0
while i<=10:
i += 1 # upper limit on wait time
T1m = a.T1
T2m = a.T2
print('T1: ' + str(a.T1) + ' T2: ' + str(a.T2))
print('Sleep 30 sec')
time.sleep(30)
if (a.T1<30 and a.T2<30 and a.T1>=T1m-0.2 and a.T2>=T2m-0.2):
break # continue when conditions met
else:
print('Not at ambient temperature')
# run step test (10 min)
for i in range(tf+1):
# set heater values
a.Q1(Q1d[i])
a.Q2(Q2d[i])
print('Time: ' + str(i) + \
' H1: ' + str(Q1d[i]) + \
' H2: ' + str(Q2d[i]) + \
' T1: ' + str(a.T1) + \
' T2: ' + str(a.T2))
# wait 1 second
time.sleep(1)
fid = open(filename,'a')
fid.write(str(i)+','+str(Q1d[i])+','+str(Q2d[i])+',' \
+str(a.T1)+','+str(a.T2)+'\n')
fid.close()
# close connection to Arduino
a.close()
# read data file
data = pd.read_csv(filename)
# plot prediction and measurement
plt.figure()
plt.subplot(2,1,1)
plt.plot(m.time,Q1.value,'r-',label='Heater 1')
plt.plot(m.time,Q2.value,'b--',label='Heater 2')
plt.ylabel('Heater (%)')
plt.legend(loc='best')
plt.subplot(2,1,2)
plt.plot(data['Time'],data['T1'],'r.',label='T1 Measured')
plt.plot(data['Time'],data['T2'],'b.',label='T2 Measured')
plt.plot(m.time,T1.value,'k-',label='T1 Predicted')
plt.plot(m.time,T2.value,'k--',label='T2 Predicted')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')
plt.xlabel('Time (sec)')
plt.savefig('tclab_dyn_pred.png')
plt.show()