如何使用数组从Gekko应用m.connection？

2024-09-27 00:20:30 发布

您现在位置：Python中文网/ 问答频道 /正文
2388 3

我试图用数组、“for cycle”和How to solve warning message in Gekko due to m.connection?的输入来解决一个机械反应模型

我不确定变量的初始化是否可以这样做。此外，我怀疑我是否正确使用了约束。先谢谢你

这是我正在编写的代码

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt 
import math as math
import pandas as pd

####measured lab data - nh2cl decay
#data set 1 
t_data1 = [0,0.08333,0.5,1,4,22.6167] #nh2cl
x_data1 = [0,4.91e-5,4.57e-5,4.74e-5,4.17e-5,2.76e-5] #8.01e-5
x_data1mgl = [0, 3.48,3.24,3.36,2.96,1.96]#5.68

#data set 2
t_data2 = [0,0.08333,0.5,1,4,22.8167]
x_data2 = [0,5.92e-5,5.7e-5,5.64e-5,5.30e-5,4.6e-5] #8.01e-5
x_data2mgl = [0,4.2,4.04,4,3.76,2.88]#5.68


####Preparing lab data and time dataframe
#combine data and time in the same dataframe
data1 = pd.DataFrame({'time':t_data1,'x1':x_data1,'x1mgl': x_data1mgl})
data2 = pd.DataFrame({'time':t_data2,'x2':x_data2,'x2mgl': x_data2mgl})

data1.set_index('time', inplace=True)
data2.set_index('time', inplace=True)

# merge dataframes
data = data1.join(data2, how='outer')

# simulation time points
dftp = pd.DataFrame({'time':np.linspace(0,5,100)})
dftp.set_index('time', inplace=True)

# merge dataframes
dff = data.join(dftp,how='outer')

z1 = (dff['x1']==dff['x1']).astype(int)
z2 = (dff['x2']==dff['x2']).astype(int)

# replace NaN with zeros
df0 = dff.fillna(value=0)


####Estimator Model
m = GEKKO(remote=False)
#m = GEKKO() 

m.time = df0.index.values

# measurements - Gekko arrays to store measured data
xm = m.Array(m.Param,2)
xm[0].value = df0['x1'].values
xm[1].value = df0['x2'].values

####inital measured values
hocl_init_val=m.Array(m.Var,2)
hocl_init_val[0].value= 8.01e-5
hocl_init_val[1].value= 8.01e-5 

nh3_init_val=m.Array(m.Var,2)
nh3_init_val[0].value= 1.37e-3 
nh3_init_val[1].value= 6.82e-4 

h_init_val=m.Array(m.Var,2)
h_init_val[0].value= 6.31e-8 
h_init_val[1].value= 2e-8

h2co3_init_val=m.Array(m.Var,2)
h2co3_init_val[0].value= 5.39e-4
h2co3_init_val[1].value= 1.88e-4

hco3_init_val=m.Array(m.Var,2)
hco3_init_val[0].value= 3.46e-3
hco3_init_val[1].value= 3.80e-3

co32_init_val=m.Array(m.Var,2)
co32_init_val[0].value= 2.16e-6
co32_init_val[1].value= 7.5e-6

alk_init_val=m.Array(m.Var,2)
alk_init_val[0].value= 3.46e-3
alk_init_val[1].value= 3.82-3

hocl=m.Array(m.Var,2)
nh3=m.Array(m.Var,2)
nh2cl=m.Array(m.Var,2,lb=1e-10)
nh2cl_meas=m.Array(m.Param,2)
nhcl2=m.Array(m.Var,2,lb=1e-10)
h=m.Array(m.Var,2)
oh=m.Array(m.Var,2)
I=m.Array(m.Var,2,lb=1e-10)
ocl=m.Array(m.Var,2)
nh4=m.Array(m.Var,2)
h2co3=m.Array(m.Var,2)
hco3=m.Array(m.Var,2)
co32=m.Array(m.Var,2)
alk=m.Array(m.Var,2)
#DOC1=m.Array(m.Var,2)
#DOC2=m.Array(m.Var,2)
cnh3=m.Array(m.Var,2)
cnh2cl=m.Array(m.Var,2)

r1=m.Array(m.Var,2)
r2=m.Array(m.Var,2)
r3=m.Array(m.Var,2)
r4=m.Array(m.Var,2)
r5=m.Array(m.Var,2)
r6=m.Array(m.Var,2)
r7=m.Array(m.Var,2)
r8=m.Array(m.Var,2)
r9=m.Array(m.Var,2)
r10=m.Array(m.Var,2)

r11=m.Array(m.Var,2)
r12=m.Array(m.Var,2)
r13=m.Array(m.Var,2)
r14=m.Array(m.Var,2)

r15=m.Array(m.Var,2)
r16=m.Array(m.Var,2)



#Define GEKKO variables that determine if time point contains data to be used in regression
#index for objective (0=not measured, 1=measured)
zm = m.Array(m.Param,2)
zm[0].value=z1
zm[1].value=z2

# fit to measurement
x=m.Array(m.Var,2)                         
x[0].value=x_data1[0]
x[1].value=x_data2[0]

meas=m.Array(m.Param,2) 

#### adjustable parameters 

kdoc1 = m.FV(1.2334e6,lb=1,ub=100000000000) #m-1h-1
kdoc2 = m.FV(3.8809e9,lb=1,ub=10000000000) #m-1h-1

x10=3.1684e-5
x20=3.1308e-5
#DOC10 = m.FV(x10,lb=1e-6,ub=1)
#DOC20 = m.FV(x20,lb=1e-6,ub=1)

#constrains
DOC10=m.Array(m.FV,2)
DOC10[0].value=x10
DOC10[1].value=x10*0.5

DOC20=m.Array(m.FV,2)
DOC20[0].value=x20
DOC20[1].value=x20*0.5

#variables connected to the parameters DOC20 and DOC10: DOC10 is DOC1 when t=0 and DOC20 is DOC2 when t=0 
DOC1=m.Array(m.SV,2,fixed_initial=False)
#DOC1[0].value=x10
#DOC1[1].value=x10*0.5

DOC2=m.Array(m.SV,2,fixed_initial=False)
#DOC2[0].value=x20
#DOC2[1].value=x20*0.5

#t = m.Param(value=m.time)

####Rate constants
k1 = m.FV(1.5e10)
k2 = m.FV(7.6e-2)
k3 = m.FV(1e6)
k4 = m.FV(2.3e-3)
k6 = m.FV(2.2e8)
k7 = m.FV(4e5)
k8 = m.FV(1e8)
k9 = m.FV(3e7)
k10 = m.FV(55)
k11 = m.FV(3.16e-8*1e10)
k12 = m.FV(1e10)
k13 = m.FV(5.01e-10*1e10)
k14 = m.FV(1e10)

k5=m.Array(m.Var,2)

for i in range(2):
    m.Connection(DOC2[i],DOC20[i],pos1=1,pos2=1,node1=1,node2=1)
    m.Connection(DOC1[i],DOC10[i],pos1=1,pos2=1,node1=1,node2=1)
    
    #variables - initial values
    hocl[i] = m.Var(value=hocl_init_val[i])
    nh3[i] = m.Var(value=nh3_init_val[i])
    nh2cl[i] = m.Var(value=x[i])
    nh2cl_meas[i] = m.Param(xm[i])
    meas[i] = m.Param(zm[i])
    nhcl2[i] = m.Var(value=0)
    h[i] = m.Var(value=h_init_val[i])
    oh[i] = m.Var(value=1e-14/h_init_val[i])
    I[i] = m.Var(value=0)
    ocl[i]= m.Var(value=0)
    nh4[i] = m.Var(value=0)
    h2co3[i] = m.Var(value=h2co3_init_val[i])
    hco3[i] = m.Var(value=hco3_init_val[i])
    co32[i] = m.Var(value=co32_init_val[i])
    alk[i] = m.Var(value=alk_init_val[i])
    DOC1[i] = m.SV(value=DOC10[i])
    DOC2[i] = m.SV(value=DOC20[i])
    cnh3[i] = m.Var(value=0)
    cnh2cl[i] = m.Var(value=0)
    

###Equations 1   
    m.Equation(k5[i] == 2.5e7*h[i]+4e4*h2co3[i]+800*hco3[i])

###Equations 2 
    m.Equation(r1[i] == k1 * hocl[i] * nh3[i])
    m.Equation(r2[i] == k2 * nh2cl[i])
    m.Equation(r3[i] == k3 * hocl[i] * nh2cl[i])
    m.Equation(r4[i] == k4 * nhcl2[i])
    m.Equation(r5[i] == k5[i] * nh2cl[i] * nh2cl[i])
    m.Equation(r6[i] == k6 * nhcl2[i] * nh3[i]* h[i])
    m.Equation(r7[i] == k7 * nhcl2[i] * oh[i])
    m.Equation(r8[i] == k8 * I[i] * nhcl2[i])
    m.Equation(r9[i] == k9 * I[i] * nh2cl[i])
    m.Equation(r10[i] == k10 * nh2cl[i] * nhcl2[i])
    
    m.Equation(r11[i] == k11*hocl[i])
    m.Equation(r12[i] == k12*h[i]*ocl[i])
    m.Equation(r13[i] == k13*nh4[i])
    m.Equation(r14[i] == k14*h[i]*nh3[i])
    
    m.Equation(r15[i] == kdoc1*DOC1[i]*nh2cl[i])
    m.Equation(r16[i] == kdoc2*DOC2[i]*hocl[i])
    
####Equations 3

    m.Equation(hocl[i].dt() == -r1[i] + r2[i] - r3[i] + r4[i] + r8[i] - r16[i] - r11[i] + r12[i])
      
    m.Equation(nh3[i].dt() == -r1[i] + r2[i] + r5[i] - r6[i] + r13[i] - r14[i])
    m.Equation(nh2cl[i].dt() == r1[i] - r2[i] - r3[i] + r4[i] - r5[i] + r6[i] - r9[i] - r10[i] - r15[i])
    m.Equation(nhcl2[i].dt() == r3[i] - r4[i] + r5[i] - r6[i] - r7[i] - r8[i] - r10[i])
    m.Equation(h[i].dt() == 0)
    
    m.Equation(oh[i] == 1e-14/h[i])
    
    m.Equation(I[i].dt() == r7[i] - r8[i] - r9[i])  
    m.Equation(ocl[i].dt() == r11[i] - r12[i])
    m.Equation(nh4[i].dt() == -r13[i]  + r14[i])
    
    m.Equation(h2co3[i] == (hco3[i]*h[i])/5.01e-7)
    
    m.Equation(hco3[i] == alk[i] - 2*co32[i] - oh[i] + h[i])
    
    m.Equation(co32[i] == (5.01e-11*hco3[i])/h[i])
    
    m.Equation(alk[i].dt() == 0)
       
    m.Equation(DOC1[i].dt() == -r15[i])
    m.Equation(DOC2[i].dt() == -r16[i])
    
    m.Equation(cnh3[i] == 17000*nh3[i])
    m.Equation(cnh2cl[i] == 70900*nh2cl[i])
    
###obj function
    
    m.Minimize(meas[i]*((nh2cl[i]-nh2cl_meas[i]))**2)
        
#Application options
m.options.SOLVER = 1 #IPOPT solver=3 #APOPT solver=1
m.options.IMODE = 8 #Dynamic Simultaneous - estimation = 5 # Dynamic sequential - estimation = 8
m.options.RTOL = 1E-3
m.options.EV_TYPE = 2 #absolute error
m.options.NODES = 2 #collocation nodes (2 up to 6)
#m.options.MAX_ITER = 500
m.open_folder()

if True:
    kdoc1.STATUS=1
    kdoc2.STATUS=1
    DOC10[i].STATUS=1
    DOC20[i].STATUS=1

m.options.TIME_SHIFT = 0
try:
    m.solve(disp=True)
except:
    print("don't stop when not finding Doc10")


print('Final SSE Objective: ' + str(m.options.objfcnval))

print('Solution')
print('kdoc1 = ' + str(kdoc1.value[0]))
print('kdoc2 = ' + str(kdoc2.value[0]))
print('DOC10_1 = ' + str(DOC10[0].value[0]))
print('DOC20_1 = ' + str(DOC20[0].value[0]))
print('DOC10_2 = ' + str(DOC10[1].value[0]))
print('DOC20_2 = ' + str(DOC20[1].value[0]))

####Graphics

plt.figure(1,figsize=(8,5))
plt.subplot(2,1,1)
plt.plot(m.time,nh2cl[0],'m',label='Predicted 1')
plt.plot(m.time,nh2cl[1],'c',label='Predicted 2')

plt.plot(t_data1,x_data1,'mo',label='Meas 1')
plt.plot(t_data2,x_data2,'co',label='Meas 2')

plt.legend(loc='best')
plt.ylabel('mol/L')
plt.legend(loc='center right', bbox_to_anchor=(1.45, 0.5), ncol=2) #shadow=True,

plt.subplot(2,1,2)
plt.plot(m.time,cnh2cl[0].value,'m',label ='C2_M1')
plt.plot(m.time,cnh2cl[1].value,'c',label ='C2_M2')

plt.plot(t_data1,x_data1mgl,'mo',label='Meas 1')
plt.plot(t_data2,x_data2mgl,'co',label='Meas 2')

plt.legend(loc='best')
plt.ylabel('mg Cl2/L')
plt.xlabel('time (h)')
plt.legend(loc='right', bbox_to_anchor=(1.4, 0.5),  ncol=2) #shadow=True,

plt.figure(2,figsize=(12,8))
plt.subplot(4,3,1)
plt.plot(m.time,hocl[i].value,label ='hocl')
plt.legend()

plt.subplot(4,3,2)
plt.plot(m.time,nh3[i].value,label ='nh3')
plt.legend()

plt.subplot(4,3,3)
plt.plot(m.time,nhcl2[i].value,label ='nhcl2')
plt.legend()

plt.subplot(4,3,4)
plt.plot(m.time,h[i].value,label ='h')
plt.legend()

plt.subplot(4,3,5)
plt.plot(m.time,oh[i].value,label ='oh')
plt.legend()

plt.subplot(4,3,6)
plt.plot(m.time,I[i].value,label ='I')
plt.legend()

plt.subplot(4,3,7)
plt.plot(m.time,ocl[i].value,label ='ocl')
plt.legend()

plt.subplot(4,3,8)
plt.plot(m.time,nh4[i].value,label ='nh4')
plt.legend()

plt.subplot(4,3,9)
plt.plot(m.time,h2co3[i].value,label ='h2co3')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,10)
plt.plot(m.time,hco3[i].value,label ='hco3')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,11)
plt.plot(m.time,co32[i].value,label ='co32')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,3,12)
plt.plot(m.time,alk[i].value,label ='alk')
plt.legend()
plt.xlabel('time (h)')

plt.show()

Tags： timeplotinitvaluevarpltvalarray
1条回答
网友
1楼 · 发布于 2024-09-27 00:20:30

18546自由度的问题维度看起来不正确。似乎只有速率常数作为自由度。您可能需要检查每个变量是否都有一个方程，并且在所有速率常数固定的情况下，它最初是以零自由度（方程数等于变量数）求解的

                                 
 APMonitor, Version 1.0.0
 APMonitor Optimization Suite
                                 
 
 
     - APM Model Size       
 Each time step contains
   Objects      :  0
   Constants    :  0
   Variables    :  141
   Intermediates:  0
   Connections  :  8
   Equations    :  36
   Residuals    :  36
 
 Variable time shift OFF
 Number of state variables:    30082
 Number of total equations: -  11536
 Number of slack variables: -  0
                    -
 Degrees of freedom       :    18546
 
                        
 Dynamic Estimation with APOPT Solver
                        
 
 Iter    Objective  Convergence
    0  1.37152E+01  4.08093E-01
    1  2.83795E+01  9.57058E-05
    2  1.36177E-01  1.62209E-05

  247  3.06822E+28  5.05213E-04
  248  3.01561E+28  5.05213E-04
  249  3.04167E+28  5.05213E-04
 
 Iter    Objective  Convergence
  250  3.06822E+28  5.05213E-04
 Maximum iterations
 
                          -
 Solver         :  APOPT (v1.0)
 Solution time  :  264.2566 sec
 Objective      :  4.856644896504421E-8
 Unsuccessful with error code  0
                          -
 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found

相关问题 更多 >

    编程相关推荐

    热门问题