Question
Combination of PSO and GEKKO: Error: Intermediate variable with no equality (=) expression
I want to optimize the best polynomial coefficients which describes a temperature profile in the interval [0, tf], where tf=100 min:
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
This temperature profile should maximize biodiesel concentration at final time (J = max x4(tf)), subject to a system of ODEs regarding to reaction velocities, and the inequality constraint: 298 K < T < 338 K.
This problem is originated from an article whose authors proposed a combination of Monte Carlo algorithm random search and Genetic Algorithm for the polynomials parameters parameterization, and Ode45 in Matlab environment for x4(tf) optimization.
Article: doi.org/10.1016/j.cherd.2021.11.001
I tried a combination of Particle Swarm Optimization and GEKKO. I want that the PSO finds the best polynomials coefficients so that GEKKO module maximizes x4(tf). So my intention is that after every PSO iteration, it computes the coefficients into the GEKKO optimization block, then the PSO algorithm verifies if that is the best objective function:
I know I have to get somewhere x4(tf) = 0.8 mol/L.
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]<G:
g=np.copy(p[i,:]) ####################
G=P[i] # registering best fobj of i
# Plotting
fig, axs = plt.subplots(2, 2, gridspec_kw={'hspace': 0.7, 'wspace': 0.7})
axs[0, 0].plot(x[:,0],x[:,1],'ro')
axs[0, 0].set_title('Iteração Inicial')
axs[0, 0].set_xlim([xmin[0], xmax[0]])
axs[0, 0].set_ylim([xmin[1], xmax[1]])
#Iterations
tmax=500
for tatual in range(tmax):
for i in range(N):
R1=np.random.uniform(0,1) # random value for R1
R2=np.random.uniform(0,1) # random value for R2
# Inertia
wmax=0.9
wmin=0.4
w=wmax-(wmax-wmin)*tatual/tmax # inertia factor
v[i,:]=w*v[i,:]+ c1*R1*(p[i,:]-x[i,:])+c2*R2*(g-x[i,:]) # velocity
x[i,:]=x[i,:]+v[i,:] # position
for d in range(D): # guarantee of bounds
if x[i,d]<xmin[d]:
x[i,d]=xmin[d]
v[i,d]=0
if x[i,d]>xmax[d]:
x[i,d]=xmax[d]
v[i,d]=0
X[i]=objective_function(x[i,:])
if X[i]<P[i]:
p[i,:]=x[i,:] # particle i best position
P[i]=X[i] # P update (best fobj)
if P[i]< G: # verify if it's better than global fobj
g=np.copy(p[i,:]) # registering best global position
G=P[i]
if tatual==49:
axs[0, 1].plot(x[:,0],x[:,1],'ro')
axs[0, 1].set_title('Iteração 20')
axs[0, 1].set_xlim([xmin[0], xmax[0]])
axs[0, 1].set_ylim([xmin[1], xmax[1]])
if tatual==99:
axs[1, 0].plot(x[:,0],x[:,1],'ro')
axs[1, 0].set_title('Iteração 100')
axs[1, 0].set_xlim([xmin[0], xmax[0]])
axs[1, 0].set_ylim([xmin[1], xmax[1]])
if tatual==499:
axs[1, 1].plot(x[:,0],x[:,1],'ro')
axs[1, 1].set_title('Iteração 499')
axs[1, 1].set_xlim([xmin[0], xmax[0]])
axs[1, 1].set_ylim([xmin[1], xmax[1]])
for ax in axs.flat:
ax.set(xlabel='x1', ylabel='x2')
print('Optimal x:', g)
print('Optimal Fobj(x):', objective_function(g))
However, I am returned with this weird GEKKO error:
Traceback (most recent call last):
exec(code, globals, locals)
P[i] = objective_function(x[i])
m.solve(disp=False)
raise Exception(apm_error)
Exception: @error: Intermediate Definition
Error: Intermediate variable with no equality (=) expression
-72.31243538-74.05913067-75.76877023-77.43020035-79.03172198
STOPPING...
I've tried both GEKKO and PSO codes separately for testing purposes, and both work. But I can't make them work together. I've tried so many things, besides I'm newby in GEKKO. If you guys could help me, I'd appreciate it a lot.
Thank you very much in advance!
