Код: Выделить всё
p0 = 0.1
p1 = (t - 50)*(3**0.5/500)
p2 = (t**2 - 100*t + 5000/3)*(3/(500000000**0.5))
T = coef0*p0 + coef1*p1 + coef2*p2
Эта проблема возникла из статьи, авторы которой предложили комбинацию алгоритма случайного поиска Монте-Карло и генетического алгоритма для параметризации параметров полиномов, а также Ode45 в среде Matlab для x4. (tf) оптимизация.
Статья: doi.org/10.1016/j.cherd.2021.11.001
Я попробовал комбинацию оптимизации Particle Swarm. и ГЕККО. Я хочу, чтобы PSO нашел лучшие коэффициенты полиномов, чтобы модуль GEKKO максимизировал x4 (tf). Поэтому я намерен, чтобы после каждой итерации PSO он вычислял коэффициенты в блоке оптимизации GEKKO, а затем алгоритм PSO проверял, является ли это лучшей целевой функцией:
Я знаю, что мне нужно получить где-то x4(tf) = 0,8 моль/л.
Код: Выделить всё
import matplotlib.pyplot as plt
import numpy as np
from gekko import GEKKO
xmin=np.array([-12,-12,-12]) # Minimum bounds
xmax=np.array([3400,3400,3400]) # Maximum bounds
# H=abs(xmax-xmin) # Diferença entre máximo e mínimo
N=30 # Particle number
# PSO parameters
c1=0.8 # individual
c2=1.2 # social
D=3 # Dimension
tmax = 500 # Maximum iterations
m = GEKKO()
def objective_function(x):
coef0, coef1, coef2 = x[0], x[1], x[2]
tf = 100 # Final time
m.time = np.linspace(0, tf, tf) # Time interval
t = m.time
# Temperature profile with coefficients to be optimized
p0 = 0.1
p1 = (t - 50)*(3**0.5/500)
p2 = (t**2 - 100*t + 5000/3)*(3/(500000000**0.5))
T = coef0*p0 + coef1*p1 + coef2*p2
# Arrhenius equations
A1, b1 = 3.92e7, 6614.83
A2, b2 = 5.77e5, 4997.98
A3, b3 = 5.88e12, 9993.96
A4, b4 = 0.98e10, 7366.64
A5, b5 = 5.35e3, 3231.18
A6, b6 = 2.15e4, 4824.87
k1 = m.Intermediate(A1*m.exp(-b1/T))
k2 = m.Intermediate(A2*m.exp(-b2/T))
k3 = m.Intermediate(A3*m.exp(-b3/T))
k4 = m.Intermediate(A4*m.exp(-b4/T))
k5 = m.Intermediate(A5*m.exp(-b5/T))
k6 = m.Intermediate(A6*m.exp(-b6/T))
# Decision variables' initial values
x1 = m.Var(value=0.3226)
x2 = m.Var(value=0)
x3 = m.Var(value=0)
x4 = m.Var(value=0)
x5 = m.Var(value=1.9356)
x6 = m.Var(value=0)
# Dynamic model
m.Equation(x1.dt() == -k1*x1*x5 + k2*x2*x4)
m.Equation(x2.dt() == k1*x1*x5 - k2*x2*x4 - k3*x2*x5 + k4*x3*x4)
m.Equation(x3.dt() == k3*x2*x5 - k4*x3*x4 - k5*x3*x5 + k6*x6*x4)
m.Equation(x4.dt() == k1*x1*x5 - k2*x2*x4 + k3*x2*x5 - k4*x3*x4 + k5*x3*x5 - k6*x6*x4)
m.Equation(x5.dt() == -(k1*x1*x5 - k2*x2*x4 + k3*x2*x5 - k4*x3*x4 + k5*x3*x5 - k6*x6*x4))
m.Equation(x6.dt() == k5*x3*x5 - k6*x6*x4)
p = np.zeros(tf)
p[-1] = 1.0
final = m.Param(value=p)
m.Maximize(x4*final)
m.options.IMODE = 6
m.solve(disp=False, debug=True)
f = x4.value[-1]
return f
#Inicialize PSO parameters
x=np.zeros((N,D))
X=np.zeros(N)
p=np.zeros((N,D)) # best position
P=np.zeros(N) # best f_obj value
v=np.zeros((N,D))
for i in range(N): # iteration for each particle
for d in range(D):
x[i,d]=xmin[d]+(xmax[d]- xmin[d])*np.random.uniform(0,1) # inicialize position
v[i,d]=0 # inicialize velocity (dx)
X[i]= objective_function(x[i,:])
p[i,:]=x[i,:]
P[i]=X[i]
if i==0:
g=np.copy(p[0,:]) ############
G=P[0] # fobj global value
if P[i]
Подробнее здесь: [url]https://stackoverflow.com/questions/78689083/combination-of-pso-and-gekko-error-intermediate-variable-with-no-equality[/url]