from gekko import GEKKO import numpy as np import matplotlib.pyplot as plt # generate training data x = np.linspace(0.0,2*np.pi,20) y = np.sin(x) # option for fitting function select = True # True / False if select: # Size with cosine function nin = 1 # inputs n1 = 1 # hidden layer 1 (linear) n2 = 1 # hidden layer 2 (nonlinear) n3 = 1 # hidden layer 3 (linear) nout = 1 # outputs else: # Size with hyperbolic tangent function nin = 1 # inputs n1 = 2 # hidden layer 1 (linear) n2 = 2 # hidden layer 2 (nonlinear) n3 = 2 # hidden layer 3 (linear) nout = 1 # outputs # Initialize gekko train = GEKKO() test = GEKKO() model = [train,test] for m in model: # input(s) m.inpt = m.Param() # layer 1 m.w1 = m.Array(m.FV, (nin,n1)) m.l1 = [m.Intermediate(m.w1[0,i]*m.inpt) for i in range(n1)] # layer 2 m.w2a = m.Array(m.FV, (n1,n2)) m.w2b = m.Array(m.FV, (n1,n2)) if select: m.l2 = [m.Intermediate(sum([m.cos(m.w2a[j,i]+m.w2b[j,i]*m.l1[j]) \ for j in range(n1)])) for i in range(n2)] else: m.l2 = [m.Intermediate(sum([m.tanh(m.w2a[j,i]+m.w2b[j,i]*m.l1[j]) \ for j in range(n1)])) for i in range(n2)] # layer 3 m.w3 = m.Array(m.FV, (n2,n3)) m.l3 = [m.Intermediate(sum([m.w3[j,i]*m.l2[j] \ for j in range(n2)])) for i in range(n3)] # output(s) m.outpt = m.CV() m.Equation(m.outpt==sum([m.l3[i] for i in range(n3)])) # flatten matrices m.w1 = m.w1.flatten() m.w2a = m.w2a.flatten() m.w2b = m.w2b.flatten() m.w3 = m.w3.flatten() # Fit parameter weights m = train m.inpt.value=x m.outpt.value=y m.outpt.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.w2a)): m.w2a[i].STATUS=1 m.w2b[i].STATUS=1 m.w2a[i].FSTATUS=1 m.w2b[i].FSTATUS=1 m.w2a[i].MEAS=1.0 m.w2b[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.SOLVER = 3 m.options.EV_TYPE = 2 m.solve(disp=False) # Test sample points m = test for i in range(len(m.w1)): m.w1[i].MEAS=train.w1[i].NEWVAL m.w1[i].FSTATUS = 1 print('w1['+str(i)+']: '+str(m.w1[i].MEAS)) for i in range(len(m.w2a)): m.w2a[i].MEAS=train.w2a[i].NEWVAL m.w2b[i].MEAS=train.w2b[i].NEWVAL m.w2a[i].FSTATUS = 1 m.w2b[i].FSTATUS = 1 print('w2a['+str(i)+']: '+str(m.w2a[i].MEAS)) print('w2b['+str(i)+']: '+str(m.w2b[i].MEAS)) for i in range(len(m.w3)): m.w3[i].MEAS=train.w3[i].NEWVAL m.w3[i].FSTATUS = 1 print('w3['+str(i)+']: '+str(m.w3[i].MEAS)) m.inpt.value=np.linspace(-2*np.pi,4*np.pi,100) m.options.IMODE = 2 m.options.SOLVER = 3 m.solve(disp=False) plt.figure() plt.plot(x,y,'bo',label='data') plt.plot(test.inpt.value,test.outpt.value,'r-',label='predict') plt.legend(loc='best') plt.ylabel('y') plt.xlabel('x') plt.show()