pythonoptimizationgekkoparticle-swarm

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!


Solution

  • Change how T is defined because it is an array of values. Also, optimization is not needed when Gekko is used for simulation.

    T = m.Param(coef0*p0 + coef1*p1 + coef2*p2,lb=298,ub=338)
    
    #m.Maximize(x4*final)
    m.options.IMODE = 4
    

    For the Temperature profile T, Gekko expects m.Param() inputs when they are vectors. Here is the new code:

    PSO Temperature Profile Optimization

    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(remote=False)
    
    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 = m.Param(coef0*p0 + coef1*p1 + coef2*p2,lb=298,ub=338)
        
        # 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 = 4
        
        m.solve(disp=False, debug=True)
        
        f = x4.value[-1]
        print(f'Objective: {f} with c0: {coef0:0.2f}, c1: {coef1:0.2f}, c2: {coef2:0.2f}')
        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))
    

    The output shows that the PSO algorithm is finding values that are within the desired range of x4(tf)>=0.8.

    Objective: 0.69040319802 with c0: 480.37, c1: 1051.17, c2: 1058.53
    Objective: 0.72953879156 with c0: 586.82, c1: 46.14, c2: 1942.71
    Objective: 0.80819933391 with c0: 2070.86, c1: 2254.40, c2: 318.32
    Objective: 0.78391848015 with c0: 597.28, c1: 178.80, c2: 3303.17
    Objective: 0.76343411822 with c0: 2237.15, c1: 505.65, c2: 959.68
    Objective: 0.76438551815 with c0: 1274.94, c1: 1470.23, c2: 1192.04
    Objective: 0.78609899052 with c0: 2079.49, c1: 910.99, c2: 1578.82
    Objective: 0.76287892509 with c0: 1204.67, c1: 1508.01, c2: 1147.28
    Objective: 0.47995105183 with c0: 243.01, c1: 865.63, c2: 136.56
    Objective: 0.79657986275 with c0: 1549.06, c1: 2994.10, c2: 1932.67
    Objective: 0.81878356134 with c0: 2615.95, c1: 2228.12, c2: 99.50
    Objective: 0.81402792602 with c0: 3301.50, c1: 3027.91, c2: 1771.54
    Objective: 0.76571967789 with c0: 1422.29, c1: 486.43, c2: 1752.17
    Objective: 0.79774007627 with c0: 1946.12, c1: 3114.17, c2: 2915.55
    Objective: 0.77935267715 with c0: 393.45, c1: 2859.27, c2: 2975.77
    Objective: 0.81096055569 with c0: 2799.38, c1: 493.08, c2: 1955.81
    Objective: 0.78843478051 with c0: 2353.52, c1: 975.26, c2: 1331.36
    Objective: 0.81830447795 with c0: 3341.08, c1: 588.18, c2: 3206.25
    Objective: 0.80886129367 with c0: 2370.03, c1: 1513.05, c2: 3226.52
    Objective: 0.70160310864 with c0: 300.87, c1: 676.89, c2: 1835.97
    Objective: 0.82196905283 with c0: 3349.32, c1: 3316.72, c2: 687.28
    Objective: 0.81945046508 with c0: 3363.60, c1: 359.99, c2: 1237.87
    Objective: 0.75740562215 with c0: 1026.78, c1: 1414.63, c2: 1488.10
    Objective: 0.81832258587 with c0: 3366.10, c1: 315.24, c2: 2954.72
    Objective: 0.75937879255 with c0: -7.19, c1: 1969.84, c2: 2659.16
    Objective: 0.76653792682 with c0: 977.31, c1: 606.03, c2: 2199.25
    Objective: 0.7927165714 with c0: 1223.03, c1: 787.54, c2: 3077.25
    ...
    

    A limitation of the PSO Temperature Profile Optimization is that the profile is limited to polynomial solutions. PSO is slow, but it does have the advantage of better finding the global optimum.

    Gekko Temperature Profile Optimization

    The other way to solve this is to use Gekko's optimization capabilities with gradient-based solvers.

    results

    import matplotlib.pyplot as plt
    import numpy as np
    from gekko import GEKKO
    
    m = GEKKO(remote=False)
    
    tf = 100 # Final time
    m.time = np.linspace(0, tf, tf) # Time interval
    t = m.time
    
    # Temperature profile
    T = m.MV(300,lb=298,ub=338)
    T.STATUS = 1
    
    # 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=True, debug=True)
    
    f = x4.value[-1]
    print(f'Objective: {f}')
    
    plt.figure(figsize=(6,4))
    plt.subplot(2,1,1)
    plt.plot(t,T,'r--',label='Temperature')
    plt.grid(); plt.legend(); plt.ylabel('T (K)')
    plt.subplot(2,1,2)
    plt.plot(t,x1,label='x1')
    plt.plot(t,x2,label='x2')
    plt.plot(t,x3,label='x3')
    plt.plot(t,x4,label='x4')
    plt.plot(t,x5,label='x5')
    plt.plot(t,x6,label='x6')
    plt.grid(); plt.legend()
    plt.xlabel('Time'); plt.ylabel('Mole frac')
    plt.tight_layout()
    plt.show()
    

    A solution of x4(tf)=0.83177 is found in 0.832 seconds.

    EXIT: Optimal Solution Found.
    
     The solution was found.
    
     The final value of the objective function is  -0.8313902000699634
    
     ---------------------------------------------------
     Solver         :  IPOPT (v3.12)
     Solution time  :  0.8324999999999999 sec
     Objective      :  -0.8313902000699634
     Successful solution
     ---------------------------------------------------
    
    Objective: 0.83177020385
    

    The solver drives the temperature to the maximum allowable temperature to maximize x4. Control the rate of change of the MV with T.DMAX=10 to limit the amount that the temperature can change each time period.

    It is also possible to optimize the values of the coefficients of a polynomial for a differentiable temperature profile. The initial point is not

    import matplotlib.pyplot as plt
    import numpy as np
    from gekko import GEKKO
    
    m = GEKKO(remote=False)
    
    tf = 100 # Final time
    m.time = np.linspace(0, tf, tf) # Time interval
    t = m.time
    
    # Temperature profile
    c0, c1, c2 = m.Array(m.FV,3,lb=-1000,ub=1e4)
    c0.STATUS = 1; c1.STATUS = 1; c2.STATUS = 1
    p0 = 0.1
    p1 = m.Param((t - 50)*(3**0.5/500))
    p2 = m.Param((t**2 - 100*t + 5000/3)*(3/(500000000**0.5)))
    T = m.Var(lb=298,ub=338)
    m.Equation(T==c0*p0 + c1*p1 + c2*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=True, debug=True)
    
    f = x4.value[-1]
    print(f'Objective: {f}')
    print(f'c0: {c0.value[0]}, c1: {c1.value[0]}, c2: {c2.value[0]}')
    
    plt.figure(figsize=(6,4))
    plt.subplot(2,1,1)
    plt.plot(t[1:],T[1:],'r--',label='Temperature')
    plt.grid(); plt.legend(); plt.ylabel('T (K)')
    plt.subplot(2,1,2)
    plt.plot(t,x1,label='x1')
    plt.plot(t,x2,label='x2')
    plt.plot(t,x3,label='x3')
    plt.plot(t,x4,label='x4')
    plt.plot(t,x5,label='x5')
    plt.plot(t,x6,label='x6')
    plt.grid(); plt.legend()
    plt.xlabel('Time'); plt.ylabel('Mole frac')
    plt.tight_layout()
    plt.savefig('results.png',dpi=300)
    plt.show()
    

    The initial condition is not calculated so it is left out of the plot.

    plt.plot(t[1:],T[1:],'r--',label='Temperature')
    

    Another interesting problem that is similar to this one is the Oil Shale Pyrolysis and other benchmark problems listed on the website.

    Oil Shale