मैं सीख रहा हूं कि पैरामीटर अनुमान समस्याओं के लिए GEKKO का उपयोग कैसे करें और पहले कदम के रूप में, मैं उदाहरण समस्याओं का विकास कर रहा हूं जो मेरे पास हैं पहले Scipy न्यूनतमकरण दिनचर्या का उपयोग करके लागू किया गया था। इनमें है APMonitor.com और में उपलब्ध जानकारी के बाद किया गया है के भीतर उपलब्ध पाठ्यक्रम। वर्तमान समस्या एक बैच रिएक्टर है एक मेथनॉल-से-हाइड्रोकार्बन प्रक्रिया का अनुकरण: http://www.daetools.com/docs/ ट्यूटोरियल्स-all.html#tutorial-che-opt-5

नीचे वर्णित कोड में मॉडल विवरण का अनुसरण किया जा सकता है, लेकिन जिन प्राथमिक चरणों पर विचार किया गया है वे हैं:

   A --> B   
   A + B --> C   
   C + B --> P   
   A --> C   
   A --> P   
   A + B --> P

जहां समय के एक फलन के रूप में ए, सी, और पी की सांद्रता के लिए प्रयोगात्मक डेटा उपलब्ध है। मॉडल का लक्ष्य छह प्राथमिक प्रतिक्रियाओं (k1-k6) के लिए दर स्थिरांक का अनुमान लगाना है। मुझे अभी जिस कठिनाई का सामना करना पड़ रहा है, वह यह है कि मेरा GEKKO मॉडल और मेरा Scipy.optimize - आधारित मॉडल एक ही प्रयोगात्मक डेटा और मापदंडों के लिए प्रारंभिक अनुमानों का उपयोग करने के बावजूद, विभिन्न पैरामीटर अनुमानों की ओर ले जाते हैं। मैंने इस मॉडल की तुलना gPROMS और एथेना विजुअल स्टूडियो का उपयोग करके विकसित किए गए मॉडल से की है, जिसमें इन बंद-सोर्स कार्यक्रमों के साथ प्राप्त पैरामीटर अनुमानों से सहमत होने वाले scipy मॉडल हैं। प्रत्येक कार्यक्रम के लिए अनुमानित पैरामीटर नीचे दिखाए गए हैं:

  • Scipy मॉडल (L-BFGS-B ऑप्टिमाइज़र): [k1 k2 k3 k4 k5 k6] = [2.779, 0., 0.197, 3.042, 2.148, 0.541]

  • GEKKO मॉडल (IPOPT अनुकूलक): [k1 k2 k3 k4 k5 k6] = [3.7766387559, 1.1826920269e-07, 0.21242442412, 4.130394645, 2.4232122905, 3.3140978171]

दिलचस्प बात यह है कि दोनों मॉडल अनुकूलन के अंत में 0.0123 के समान उद्देश्य फ़ंक्शन मान की ओर ले जाते हैं और प्रजातियों की एकाग्रता बनाम समय के भूखंडों में समान दिखते हैं। मैंने GEKKO के ऑप्टिमाइज़र को बदलने की कोशिश की है और बिना किसी लाभ के 1E-8 तक सहनशीलता को कड़ा किया है। मेरा अनुमान है कि मेरा GEKKO मॉडल ठीक से स्थापित नहीं है, लेकिन मुझे इसके साथ समस्या नहीं मिल रही है। मॉडल विसंगतियों के लिए अग्रणी संभावित मुद्दों पर मुझे इंगित करने में किसी भी मदद की सराहना की जाएगी। मैं नीचे दो लिपियों को संलग्न करता हूं:

स्किपी मॉडल

import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import minimize
import matplotlib.pyplot as plt

#Experimental data
times  = np.array([0.0, 0.071875, 0.143750, 0.215625, 0.287500, 0.359375, 0.431250,
                      0.503125, 0.575000, 0.646875, 0.718750, 0.790625, 0.862500,
                      0.934375, 1.006250, 1.078125, 1.150000])
A_obs = np.array([1.0, 0.552208, 0.300598, 0.196879, 0.101175, 0.065684, 0.045096,
                      0.028880, 0.018433, 0.011509, 0.006215, 0.004278, 0.002698,
                      0.001944, 0.001116, 0.000732, 0.000426])
C_obs = np.array([0.0, 0.187768, 0.262406, 0.350412, 0.325110, 0.367181, 0.348264,
                      0.325085, 0.355673, 0.361805, 0.363117, 0.327266, 0.330211,
                      0.385798, 0.358132, 0.380497, 0.383051])
P_obs = np.array([0.0, 0.117684, 0.175074, 0.236679, 0.234442, 0.270303, 0.272637,
                      0.274075, 0.278981, 0.297151, 0.297797, 0.298722, 0.326645,
                      0.303198, 0.277822, 0.284194, 0.301471])

def rxn(x, k): #rate equations in power law form r = k [A][B]
    A = x[0]
    B = x[1]
    C = x[2]
    P = x[3]
    
    k1 = k[0]
    k2 = k[1]
    k3 = k[2]
    k4 = k[3]
    k5 = k[4]
    k6 = k[5]
    
    r1 = k1 * A
    r2 = k2 * A * B
    r3 = k3 * C * B
    r4 = k4 * A
    r5 = k5 * A
    r6 = k6 * A * B
    
    return [r1, r2, r3, r4, r5, r6] #returns reaction rate of each equation

#mass balance diff eqs, function calls rxn function 

def mass_balances(t, x, *args): 
        k = args
        r = rxn(x, k)
        dAdt = - r[0] - r[1] - r[3] - r[4] - r[5]
        dBdt = + r[0] - r[1] - r[2] - r[5]
        dCdt = + r[1] - r[2] + r[3]
        dPdt = + r[2] + r[4] + r[5]

        return [dAdt, dBdt, dCdt, dPdt]
    
IC = [1.0, 0, 0, 0] #Initial conditions of species A, B, C, P
ki= [1, 1, 1, 1, 1, 1]

#Objective function definition

def obj_fun(k):   
#solve initial value problem over time span of data
    sol  = solve_ivp(mass_balances,[min(times),max(times)],IC, args = (k), t_eval=(times)) 
    y_model = np.vstack((sol.y[0],sol.y[2],sol.y[3])).T
    obs = np.vstack((A_obs, C_obs, P_obs)).T
    err = np.sum((y_model-obs)**2)
   
    return err

bnds = ((0, None), (0, None),(0, None),(0, None),(0, None),(0, None))
model = minimize(obj_fun,ki, bounds=bnds, method = 'L-BFGS-B')
k_opt = model.x

print(k_opt.round(decimals = 3))

y_calc = solve_ivp(mass_balances,[min(times),max(times)],IC, args = (model.x), t_eval=(times)) 

plt.plot(y_calc.t, y_calc.y.T)
plt.plot(times,A_obs,'bo')
plt.plot(times,C_obs,'gx')
plt.plot(times,P_obs,'rs')

GEKKO मॉडल

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

#Experimental data
times  = np.array([0.0, 0.071875, 0.143750, 0.215625, 0.287500, 0.359375, 0.431250,
                      0.503125, 0.575000, 0.646875, 0.718750, 0.790625, 0.862500,
                      0.934375, 1.006250, 1.078125, 1.150000])
A_obs = np.array([1.0, 0.552208, 0.300598, 0.196879, 0.101175, 0.065684, 0.045096,
                      0.028880, 0.018433, 0.011509, 0.006215, 0.004278, 0.002698,
                      0.001944, 0.001116, 0.000732, 0.000426])
C_obs = np.array([0.0, 0.187768, 0.262406, 0.350412, 0.325110, 0.367181, 0.348264,
                      0.325085, 0.355673, 0.361805, 0.363117, 0.327266, 0.330211,
                      0.385798, 0.358132, 0.380497, 0.383051])
P_obs = np.array([0.0, 0.117684, 0.175074, 0.236679, 0.234442, 0.270303, 0.272637,
                      0.274075, 0.278981, 0.297151, 0.297797, 0.298722, 0.326645,
                      0.303198, 0.277822, 0.284194, 0.301471])


m = GEKKO(remote = False)

t = m.time = times


Am = m.CV(value=A_obs, lb = 0)
Cm = m.CV(value=C_obs, lb = 0)
Pm = m.CV(value=P_obs, lb = 0)

A = m.Var(1, lb = 0)
B = m.Var(0, lb = 0)
C = m.Var(0, lb = 0)
P = m.Var(0, lb = 0)

Am.FSTATUS = 1
Cm.FSTATUS = 1
Pm.FSTATUS = 1
    
k1 = m.FV(1, lb = 0)
k2 = m.FV(1, lb = 0)
k3 = m.FV(1, lb = 0)
k4 = m.FV(1, lb = 0)
k5 = m.FV(1, lb = 0)
k6 = m.FV(1, lb = 0)

k1.STATUS = 1
k2.STATUS = 1
k3.STATUS = 1
k4.STATUS = 1
k5.STATUS = 1
k6.STATUS = 1

r1 = m.Var(0, lb = 0)
r2 = m.Var(0, lb = 0)
r3 = m.Var(0, lb = 0)
r4 = m.Var(0, lb = 0)
r5 = m.Var(0, lb = 0)
r6 = m.Var(0, lb = 0)
   
m.Equation(r1 == k1 * A)
m.Equation(r2 == k2 * A * B)
m.Equation(r3 == k3 * C * B)
m.Equation(r4 == k4 * A)
m.Equation(r5 == k5 * A)
m.Equation(r6 == k6 * A * B)
    

#mass balance diff eqs, function calls rxn function 
m.Equation(A.dt() == - r1 - r2 - r4 - r5 - r6)
m.Equation(B.dt() ==  r1 - r2 - r3 - r6)
m.Equation(C.dt() ==  r2 - r3 + r4)
m.Equation(P.dt() ==  r3 + r5 + r6)

m.Obj((A-Am)**2+(P-Pm)**2+(C-Cm)**2)


m.options.IMODE = 5
m.options.SOLVER = 3 #IPOPT optimizer
m.options.RTOL = 1E-8
m.options.OTOL = 1E-8
m.solve()

k_opt = [k1.value[0],k2.value[0], k3.value[0], k4.value[0], k5.value[0], k6.value[0]]
print(k_opt)
plt.plot(t,A)
plt.plot(t,C)
plt.plot(t,P)
plt.plot(t,B)
plt.plot(times,A_obs,'bo')
plt.plot(times,C_obs,'gx')
plt.plot(times,P_obs,'rs')
2
juanmvenegas 11 अगस्त 2020, 03:57

1 उत्तर

सबसे बढ़िया उत्तर

यहां कुछ सुझाव दिए गए हैं:

  • बेहतर एकीकरण सटीकता प्राप्त करने के लिए m.options.NODES=3 या इससे अधिक 6 तक सेट करें।
  • चर के बजाय पैरामीटर के रूप में Am, Cm, Pm सेट करें। वे निश्चित इनपुट हैं।
  • विभिन्न प्रारंभिक स्थितियों का प्रयास करें। कई स्थानीय मिनीमा हो सकते हैं।
  • उद्देश्य फलन समतल हो सकता है ताकि विभिन्न पैरामीटर मान समान उद्देश्य फलन मान दें। आप पैरामीटर कॉन्फिडेंस इंटरवल की जांच कर सकते हैं यह देखने के लिए कि डेटा संकीर्ण या विस्तृत संयुक्त विश्वास क्षेत्र।

यहां संशोधनों के साथ परिणाम दिए गए हैं:

parameter estimation

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

#Experimental data
times  = np.array([0.0, 0.071875, 0.143750, 0.215625, 0.287500, 0.359375, 0.431250,
                      0.503125, 0.575000, 0.646875, 0.718750, 0.790625, 0.862500,
                      0.934375, 1.006250, 1.078125, 1.150000])
A_obs = np.array([1.0, 0.552208, 0.300598, 0.196879, 0.101175, 0.065684, 0.045096,
                      0.028880, 0.018433, 0.011509, 0.006215, 0.004278, 0.002698,
                      0.001944, 0.001116, 0.000732, 0.000426])
C_obs = np.array([0.0, 0.187768, 0.262406, 0.350412, 0.325110, 0.367181, 0.348264,
                      0.325085, 0.355673, 0.361805, 0.363117, 0.327266, 0.330211,
                      0.385798, 0.358132, 0.380497, 0.383051])
P_obs = np.array([0.0, 0.117684, 0.175074, 0.236679, 0.234442, 0.270303, 0.272637,
                      0.274075, 0.278981, 0.297151, 0.297797, 0.298722, 0.326645,
                      0.303198, 0.277822, 0.284194, 0.301471])

m = GEKKO(remote=False)

t = m.time = times

Am = m.Param(value=A_obs, lb = 0)
Cm = m.Param(value=C_obs, lb = 0)
Pm = m.Param(value=P_obs, lb = 0)

A = m.Var(1, lb = 0)
B = m.Var(0, lb = 0)
C = m.Var(0, lb = 0)
P = m.Var(0, lb = 0)

k = m.Array(m.FV,6,value=1,lb=0)  
for ki in k:
    ki.STATUS = 1
k1,k2,k3,k4,k5,k6 = k

r1 = m.Var(0, lb = 0)
r2 = m.Var(0, lb = 0)
r3 = m.Var(0, lb = 0)
r4 = m.Var(0, lb = 0)
r5 = m.Var(0, lb = 0)
r6 = m.Var(0, lb = 0)
   
m.Equation(r1 == k1 * A)
m.Equation(r2 == k2 * A * B)
m.Equation(r3 == k3 * C * B)
m.Equation(r4 == k4 * A)
m.Equation(r5 == k5 * A)
m.Equation(r6 == k6 * A * B)

#mass balance diff eqs, function calls rxn function 
m.Equation(A.dt() == - r1 - r2 - r4 - r5 - r6)
m.Equation(B.dt() ==  r1 - r2 - r3 - r6)
m.Equation(C.dt() ==  r2 - r3 + r4)
m.Equation(P.dt() ==  r3 + r5 + r6)

m.Minimize((A-Am)**2)
m.Minimize((P-Pm)**2)
m.Minimize((C-Cm)**2)

m.options.IMODE = 5
m.options.SOLVER = 3 #IPOPT optimizer
m.options.RTOL = 1E-8
m.options.OTOL = 1E-8
m.options.NODES = 5
m.solve()

k_opt = []
for ki in k:
    k_opt.append(ki.value[0])
print(k_opt)

plt.plot(t,A)
plt.plot(t,C)
plt.plot(t,P)
plt.plot(t,B)
plt.plot(times,A_obs,'bo')
plt.plot(times,C_obs,'gx')
plt.plot(times,P_obs,'rs')
plt.show()
2
John Hedengren 11 अगस्त 2020, 03:30