How to smooth a MV function in GEKKO?

I'm working on a dynamic optimization problem very similar with this one: Model Reduction with Intermediates in GEKKO. I want to find the best temperature profile that maximizes Biodiesel concentration at final time (x4(tf)), and to compare to a polynomial temperature profile.

I would like to make it smooth and differentiable. Utilizing only T = m.MV(value=298,ub=398,lb=298) and T.STATUS = 1, it's not realistic for a state variable. I tried playing with all sorts of values for T.DCOST, T.DMAX, MV_STEP_HOR, ... but got to nowhere, not even similar to the example above. Is there another tuning parameter I'm missing here?

Best regards, and thank you very much in advance!

Here is my code:

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

m = GEKKO()

tf = 100 # Tempo final # Final time
m.time = np.linspace(0, tf, tf) # Intervalo de tempo # Time interval

# Temperatura
T = m.MV(value=298, ub=338, lb=298)
T.STATUS = 1
T.DCOST = 1e-5
T.DMAX = 0.7

#print(T.DCOST)

# Variáveis de decisão
# Decision variables
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)


# Constantes da Equação de Arrhenius para velocidade de reação
# Arrhenius Equation
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))


# Modelo dinâmico
# 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)

print(f'x4(tf): {x4.value[-1]:0.4f}')
print(T.value[-1])

# Plot de concentração e perfil de temperatura
# Concentration and temperature plot
plt.figure(figsize=(7,4))
plt.subplot(2,1,1)
plt.plot(m.time, x4.value, label='Biodiesel')
plt.xlabel('Tempo (min)')
plt.ylabel('Concentração (mol/L)')
plt.legend(loc='best')
plt.grid(True)

plt.subplot(2,1,2)
plt.plot(m.time, T.value, label=f'GEKKO DCOST={T.DCOST}', c='r')
# plt.axhline(y=338, c='y', linestyle='--')
# plt.axhline(y=298, c='y', linestyle='--')
plt.xlabel('Tempo (min)')
plt.ylabel('Temperatura (K)')
plt.legend(loc='best')
plt.grid(True); plt.tight_layout()
plt.savefig('biodiesel.png',dpi=300)
plt.show()

The MV variable is adjustable at every time point in the horizon with either a Zero-Order Hold (Steps) or First-Order Hold (Ramp). In this case, the MV saturates at the upper temperature limit because of the upper bound of 338.

T = m.MV(value=298, ub=338, lb=298)

To make it differentiable, the MV variable would need to become a polynomial function or other function that is differentiable. Coefficients from the function can be adjusted to create the MV profile that optimizes the biodiesel production.

# Temperatura
T0 = 298
T = m.Var(value=T0,lb=298,ub=338)
tm = m.Param(m.time)
c1,c2,c3 = m.Array(m.FV,3,lb=-1e5,ub=1e5)
c1.STATUS=1; c2.STATUS=1; c3.STATUS=1
m.Equation(T==T0 + c1*tm + c2*tm**2 + c3*tm**3)

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

m = GEKKO()

tf = 100 # Tempo final # Final time
m.time = np.linspace(0, tf, tf) # Intervalo de tempo # Time interval

# Temperatura
T0 = 298
T = m.Var(value=T0,lb=298,ub=338)
tm = m.Param(m.time)
c1,c2,c3 = m.Array(m.FV,3,lb=-1e5,ub=1e5)
c1.STATUS=1; c2.STATUS=1; c3.STATUS=1
m.Equation(T==T0 + c1*tm + c2*tm**2 + c3*tm**3)

# Variáveis de decisão
# Decision variables
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)

# Constantes da Equação de Arrhenius para velocidade de reação
# Arrhenius Equation
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))

# Modelo dinâmico
# 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)

print(f'x4(tf): {x4.value[-1]:0.4f}')
print(T.value[-1])

# Plot de concentração e perfil de temperatura
# Concentration and temperature plot
plt.figure(figsize=(7,4))
plt.subplot(2,1,1)
plt.plot(m.time, x4.value, label='Biodiesel')
plt.xlabel('Tempo (min)')
plt.ylabel('Concentração (mol/L)')
plt.legend(loc='best')
plt.grid(True)

plt.subplot(2,1,2)
plt.plot(m.time, T.value, 'r--', label='Temperatura')
plt.xlabel('Tempo (min)')
plt.ylabel('Temperatura (K)')
plt.legend(loc='best')
plt.grid(True); plt.tight_layout()
plt.savefig('biodiesel.png',dpi=300)
plt.show()

However, the constraint of a differentiable function may decrease the overall objective because of the fewer degrees of freedom (3 decision variables of the 3rd order polynomial versus 99 decision variables that are adjusted for every time step with the MV type).