Оценка параметров устойчивого состояния для CSTR с использованием GEKKO - PullRequest
Новая
       9

Оценка параметров устойчивого состояния для CSTR с использованием GEKKO

3 голосов
/ 02 ноября 2019

Я бы хотел подогнать константы реакции (k0 и EoverR), используя данные о концентрации CSTR в стационарном состоянии относительно температуры реактора. Я использовал простое линейное уравнение, чтобы сгенерировать данные операции, чтобы соответствовать. (Ca_data = -1,5 * T_reactor / 100 + 4,2)

Поскольку это стационарные данные, не требуется никакого временного шага (m.time). Пожалуйста, дайте совет, как преобразовать приведенный ниже код моделирования в код оценки для «Ca vs. T_reactor». 

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

# Feed Temperature (K)
Tf = 350
# Feed Concentration (mol/m^3)
Caf = 1

# Steady State Initial Conditions for the States
Ca_ss = 1
T_ss = 304

#%% GEKKO
m = GEKKO(remote=True)
m.time = np.linspace(0, 25, 251)

# Volumetric Flowrate (m^3/sec)
q = 100
# Volume of CSTR (m^3)
V = 100
# Density of A-B Mixture (kg/m^3)
rho = 1000
# Heat capacity of A-B Mixture (J/kg-K)
Cp = 0.239
# Heat of reaction for A->B (J/mol)
mdelH = 5e4
# E - Activation energy in the Arrhenius Equation (J/mol)
# R - Universal Gas Constant = 8.31451 J/mol-K
EoverR = 8700
# Pre-exponential factor (1/sec)
k0 = 3.2e15
# U - Overall Heat Transfer Coefficient (W/m^2-K)
# A - Area - this value is specific for the U calculation (m^2)
UA = 5e4

# initial conditions = 280
T0 = 304
Ca0 = 1.0

T = m.MV(value=T_ss)
rA = m.Var(value=0)
Ca = m.CV(value=Ca_ss)

m.Equation(rA == k0*m.exp(-EoverR/T)*Ca)
m.Equation(Ca.dt() == q/V*(Caf - Ca) - rA)

m.options.IMODE = 1
m.options.SOLVER = 3

T_reactor = np.linspace(220, 260, 11)
Ca_results = np.zeros(np.size(T_reactor))
for i in range(np.size(T_reactor)):
    T.Value = T_reactor[i]
    m.solve(disp=True)
    Ca_results[i] = Ca[-1]

Ca_data = -1.5*T_reactor/100 + 4.2 # for generating the operation data

# Plot the results
plt.plot(T_reactor,Ca_data,'bo',linewidth=3)
plt.plot(T_reactor,Ca_results,'r-',linewidth=3)
plt.ylabel('Ca (mol/L)')
plt.xlabel('Temperature (K)')
plt.legend(['Reactor Concentration'],loc='best')
plt.show()

1 Ответ

2 голосов
/ 04 ноября 2019

В Gekko существует режим оценки установившегося состояния (IMODE=2) для линейной или нелинейной регрессии. Двумя примерами являются нелинейная регрессия и регрессия цены на энергию . Для опубликованной проблемы, вот несколько рекомендаций:

  • Решите проблему регрессии одним решением, а не в цикле. Таким образом, выбранные вами параметры будут соответствовать всему диапазону, а не только одной точке.
  • Определите параметры, которые необходимо отрегулировать, чтобы минимизировать ошибку между данными и предсказаниями модели. Это должен быть тип m.FV() для параметров, которые имеют одно значение, и m.MV() для параметров, которые имеют разные значения для каждой временной точки.
  • Установите Ca.FSTATUS=1, чтобы сообщить решателю, что он должен попытаться сопоставитьпредсказания Ca для данных, загруженных в Ca.value.
  • Установите kf.STATUS=1, чтобы сообщить решающему, что это параметр, который должен быть скорректирован для минимизации ошибки Ca.
  • Необязательно: установите kf настраиваемый параметр вместо k0 напрямую, чтобы улучшить масштабирование проблемы. Большие значения (например,> 1e10 или <-1e10) могут создавать проблемы для решателя (без автоматического масштабирования), поскольку градиенты малы. Я создал новый параметр <code>kf в качестве промежуточного звена, чтобы также удалить дополнительное уравнение.

Regression Fit

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO(remote=True)
Tf = 350
Caf = 1
q = 100
V = 100
rho = 1000
Cp = 0.239
mdelH = 5e4
EoverR = 8700
k0 = 3.2e15
UA = 5e4
T = m.MV()
Ca = m.CV()

# new parameter to estimate
kf = m.FV(1,lb=0.5,ub=2.0)
kf.STATUS = 1

rA = m.Intermediate(kf*k0*m.exp(-EoverR/T)*Ca)
m.Equation(Ca.dt() == q/V*(Caf - Ca) - rA)

m.options.IMODE = 2
m.options.SOLVER = 3

# generate data
T_reactor = np.linspace(220, 260, 11)
Ca_data = -1.5*T_reactor/100 + 4.2

# insert data
T.value = T_reactor
Ca.value = Ca_data
Ca.FSTATUS = 1 # fit Ca

m.solve()

print('kf = ' + str(kf.value[0]))
print('k = ' + str(kf.value[0]*k0))

# Plot the results
plt.plot(T_reactor,Ca_data,'bo',linewidth=3)
plt.plot(T_reactor,Ca.value,'r-',linewidth=3)
plt.ylabel('Ca (mol/L)')
plt.xlabel('Temperature (K)')
plt.legend(['Reactor Concentration','Regression Fit'],loc='best')
plt.show()

Вы можете выбрать любое количество параметров для оценки, чтобы улучшить подгонку. Это не должно быть ограничено только kf. В вашем сообщении упоминается, что EoverR является еще одним потенциальным параметром для оценки, но это может существенно не улучшить подгонку, поскольку k0 и EoverR являются ко-линейными. Оба параметра могут быть увеличены или уменьшены и дают почти одинаковое решение. В качестве предостережения, необходимо значительное изменение температуры, чтобы оценить оба.

...