如何解决由于 m.connection,如何解决 Gekko 中的警告消息?
我正在使用 m.connection 来估计变量的初始条件,但我收到了 12 条警告消息,例如:
我不知道如何解决这些消息。
我遵循这个解释“如果 pos1 或 pos2 不是 None,则关联的 var 必须是 GEKKO 变量,并且位置是变量的(0 索引)时间离散索引”来写 m.Connection(var1,var2,pos1=None,pos2=None,node1='end',node2='end')。 https://gekko.readthedocs.io/en/latest/quick_start.html#connections
提前致谢。
$ShardMapManager = Get-ShardMapManager -Username $DbUsrId -Password $DbPwd
解决方法
底层节点结构有一个 1-index 而不是 Python 中常见的 0-index。使用 pos1=1
和 pos2=1
解决警告。
m.Connection(cnf1,cnf01,pos1=1,pos2=1,node1=1,node2=1)
m.Connection(cnf2,cnf02,node2=1)
m.Connection(cnf3,cnf03,node2=1)
m.Connection(cnf4,cnf04,node2=1)
m.Connection(cnf5,cnf05,node2=1)
m.Connection(cnf6,cnf06,node2=1)
另一个问题是 Gekko 变量通常不应该用于初始化其他值。我建议设置 x0=1.3
并使用该浮点数来初始化变量。将 m.Var()
更改为 m.SV()
以避免在连接期间将 m.Var()
重新分类为 m.FV()
。 m.SV()
是一种提升类型的变量,与 m.FV()
具有相同的优先级。这是一个完整的脚本,尽管结果看起来并不理想。
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import math as math
import pandas as pd
tm1 = [0,0.0667,0.5,1,4,22.61667]
mca1 = [5.68,3.48,3.24,3.36,2.96,1.96]
tm2 = [0,0.08333,4.25,22.8167]
mca2 = [5.68,4.20,4.04,4.00,3.76,2.88]
tm3 = [0,4.33,22.9500]
mca3 = [5.68,4.64,4.52,4.56,4.24,3.72]
tm4 = [0,4.0833,23.0833]
mca4 =[18.90,15.4,14.3,15.10,13.50,10.90]
tm5 = [0,4.5,23.2167]
mca5 =[18.90,15.5,16.30,16,14.70,13.00]
tm6 = [0,4.6667,23.3333 ]
mca6 = [18.90,15.8,11.70,12,9.5 ]
df1=pd.DataFrame({'time':tm1,'x1':mca1})
df2=pd.DataFrame({'time':tm2,'x2':mca2})
df3=pd.DataFrame({'time':tm3,'x3':mca3})
df4=pd.DataFrame({'time':tm4,'x4':mca4})
df5=pd.DataFrame({'time':tm5,'x5':mca5})
df6=pd.DataFrame({'time':tm6,'x6':mca6})
df1.set_index('time',inplace=True)
df2.set_index('time',inplace=True)
df3.set_index('time',inplace=True)
df4.set_index('time',inplace=True)
df5.set_index('time',inplace=True)
df6.set_index('time',inplace=True)
#simulation time points
dfx = pd.DataFrame({'time':np.linspace(0,25,101)})
dfx.set_index('time',inplace=True)
#merge dataframes
dfxx=dfx.join(df1,how='outer')
dfxxx=dfxx.join(df2,how='outer')
dfxxxx=dfxxx.join(df3,how='outer')
dfxxxxx=dfxxxx.join(df4,how='outer')
dfxxxxxx=dfxxxxx.join(df5,how='outer')
df=dfxxxxxx.join(df6,how='outer')
# get True (1) or False (0) for measurement
df['meas1']=(df['x1'].values==df['x1'].values).astype(int)
df['meas2']=(df['x2'].values==df['x2'].values).astype(int)
df['meas3']=(df['x3'].values==df['x3'].values).astype(int)
df['meas4']=(df['x4'].values==df['x4'].values).astype(int)
df['meas5']=(df['x5'].values==df['x5'].values).astype(int)
df['meas6']=(df['x6'].values==df['x6'].values).astype(int)
#replace NaN with zeros
df0=df.fillna(value=0)
m = GEKKO()
m.time = df0.index.values
meas1 = m.Param(df0['meas1'].values)
meas2 = m.Param(df0['meas2'].values)
meas3 = m.Param(df0['meas3'].values)
meas4 = m.Param(df0['meas4'].values)
meas5 = m.Param(df0['meas5'].values)
meas6 = m.Param(df0['meas6'].values)
#adjustable Parameters
kf=m.FV(1.3,lb=0.01,ub=10)
ks=m.FV(1.3,ub=10)
x0 = 1.3
cnf01=m.FV(x0,ub=10)
cns01=m.FV(x0,ub=10)
#constrains
cnf02=m.FV(value=x0*0.5,lb=x0*0.5,ub=x0*0.5)
cns02=m.FV(value=x0*0.5,ub=x0*0.5)
cnf03=m.FV(value=x0*0.25,lb=x0*0.25,ub=x0*0.25)
cns03=m.FV(value=x0*0.25,ub=x0*0.25)
cnf04=m.FV(value=x0,lb=x0,ub=x0)
cns04=m.FV(value=x0,ub=x0)
cnf05=m.FV(value=x0*0.5,ub=x0*0.5)
cns05=m.FV(value=x0*0.5,ub=x0*0.5)
cnf06=m.FV(value=x0*0.25,ub=x0*0.25)
cns06=m.FV(value=x0*0.25,ub=x0*0.25)
#Variables
c1 = m.SV(value=mca1[0])
c2 = m.SV(value=mca2[0])
c3 = m.SV(value=mca3[0])
c4 = m.SV(value=mca4[0])
c5 = m.SV(value=mca5[0])
c6 = m.SV(value=mca6[0])
cm1 = m.Param(df0['x1'].values)
cm2 = m.Param(df0['x2'].values)
cm3 = m.Param(df0['x3'].values)
cm4 = m.Param(df0['x4'].values)
cm5 = m.Param(df0['x5'].values)
cm6 = m.Param(df0['x6'].values)
m.Minimize((meas1*(c1-cm1)**2)+(meas2*(c2-cm2)**2)\
+(meas3*(c3-cm3)**2)+(meas4*(c4-cm4)**2)\
+(meas5*(c5-cm5)**2)+(meas6*(c6-cm6)**2))
cnf1=m.SV(value=x0,fixed_initial=False)
cns1=m.SV(value=x0,fixed_initial=False)
cnf2=m.SV(value=x0,fixed_initial=False)
cns2=m.SV(value=x0,fixed_initial=False)
cnf3=m.SV(value=x0,fixed_initial=False)
cns3=m.SV(value=x0,fixed_initial=False)
cnf4=m.SV(value=x0,fixed_initial=False)
cns4=m.SV(value=x0,fixed_initial=False)
cnf5=m.SV(value=x0,fixed_initial=False)
cns5=m.SV(value=x0,fixed_initial=False)
cnf6=m.SV(value=x0,fixed_initial=False)
cns6=m.SV(value=x0,fixed_initial=False)
#Equations
t = m.Param(value=m.time)
m.Connection(cnf1,node2=1)
m.Connection(cns1,cns01,node2=1)
m.Connection(cns2,cns02,node2=1)
m.Connection(cns3,cns03,node2=1)
m.Connection(cns4,cns04,node2=1)
m.Connection(cns5,cns05,node2=1)
m.Connection(cns6,cns06,node2=1)
m.Equation(cnf1.dt()==-kf*c1*cnf1)
m.Equation(cns1.dt()==-ks*c1*cns1)
m.Equation(c1.dt()==cnf1.dt()+cns1.dt())
m.Equation(cnf2.dt()==-kf*c2*cnf2)
m.Equation(cns2.dt()==-ks*c2*cns2)
m.Equation(c2.dt()==cnf2.dt()+cns2.dt())
m.Equation(cnf3.dt()==-kf*c3*cnf3)
m.Equation(cns3.dt()==-ks*c3*cns3)
m.Equation(c3.dt()==cnf3.dt()+cns3.dt())
m.Equation(cnf4.dt()==-kf*c4*cnf4)
m.Equation(cns4.dt()==-ks*c4*cns4)
m.Equation(c4.dt()==cnf4.dt()+cns4.dt())
m.Equation(cnf5.dt()==-kf*c5*cnf5)
m.Equation(cns5.dt()==-ks*c5*cns5)
m.Equation(c5.dt()==cnf5.dt()+cns5.dt())
m.Equation(cnf6.dt()==-kf*c6*cnf6)
m.Equation(cns6.dt()==-ks*c6*cns6)
m.Equation(c6.dt()==cnf6.dt()+cns6.dt())
#Options
m.options.SOLVER = 1 # APOPT solver
m.options.IMODE = 5 # Dynamic Simultaneous - estimation = MHE
m.options.EV_TYPE = 2 # Squared error
m.options.NODES = 3 # Collocation nodes (2,5)
if True:
kf.STATUS=1
ks.STATUS=1
cnf01.STATUS=1
cns01.STATUS=1
cnf02.STATUS=1
cns02.STATUS=1
cnf03.STATUS=1
cns03.STATUS=1
cnf04.STATUS=1
cns04.STATUS=1
cnf05.STATUS=1
cns05.STATUS=1
cnf06.STATUS=1
cns06.STATUS=1
m.options.TIME_SHIFT = 0
try:
m.solve(disp=True)
except:
print("don't stop when not finding cnf01...cnf06")
#m.open_folder()
print('Final SSE Objective: ' + str(m.options.objfcnval))
print('Solution')
print('cnf01 = ' + str(cnf1.value[0]))
print('cns01 = ' + str(cns1.value[0]))
print('kf = ' + str(kf.value[0]))
print('ks = ' + str(ks.value[0]))
print('cns02 = '+ str(cns2.value[0]))
print('cnf02 = '+ str(cnf2.value[0]))
print('cns03 = '+ str(cns3.value[0]))
print('cnf03 = '+ str(cnf3.value[0]))
print('cns04 = '+ str(cns4.value[0]))
print('cnf04 = '+ str(cnf4.value[0]))
print('cns05 = '+ str(cns5.value[0]))
print('cnf05 = '+ str(cnf5.value[0]))
print('cns06 = '+ str(cns6.value[0]))
print('cnf06 = '+ str(cnf6.value[0]))
plt.figure(1,figsize=(8,5))
plt.plot(m.time,c1.value,'r',label='Predicted c1')
plt.plot(m.time,c2.value,'y',label='Predicted c2')
plt.plot(m.time,c3.value,'c',label='Predicted c3')
plt.plot(m.time,c4.value,'g',label='Predicted c4')
plt.plot(m.time,c5.value,'b',label='Predicted c5')
plt.plot(m.time,c6.value,'m',label='Predicted c6')
plt.plot(tm1,mca1,'rx',label='Meas c1')
plt.plot(tm2,mca2,'yx',label='Meas c2')
plt.plot(tm3,mca3,'cx',label='Meas c3')
plt.plot(tm4,mca4,'go',label='Meas c4')
plt.plot(tm5,mca5,'bo',label='Meas c5')
plt.plot(tm6,mca6,'mo',label='Meas c6')
plt.xlabel('time (h)')
plt.ylabel('Concentration (mgCl2/L)')
plt.legend(loc='upper center',bbox_to_anchor=(0.5,-0.15),ncol=2)
plt.show()
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。