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()