微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

由于 m.connection,如何解决 Gekko 中的警告消息?

如何解决由于 m.connection,如何解决 Gekko 中的警告消息?

我正在使用 m.connection 来估计变量的初始条件,但我收到了 12 条警告消息,例如:

enter image description here

此外,APM 文件显示

enter image description here

我不知道如何解决这些消息。

我遵循这个解释“如果 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=1pos2=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() 具有相同的优先级。这是一个完整的脚本,尽管结果看起来并不理想。

results

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 举报,一经查实,本站将立刻删除。