Defining previous-timestep dependent ramping constraints in Python GEKKO

I am building an MPC optimization (IMODE = 6) in Python GEKKO with multiple state (SV) and manipulated (MV) variables. Inside my problem, I would like to define ramping constraints of the following kind:

|x(tk) - x(tk-1)| <= Cx for the SV, where Cx is a constant, tk is the current time step and tk-1 is the previous

0.9 <= u(tk) / u(tk-1) <= 1.1 for the MV

I tried using the build-in tuning parameter DCOST for the MV, but as far as I understand it suggests using a constant value, which is not my case. Alternatively I tried solving the model iteratively with 1 time step per iteration and taking the values of the previous iteration as starting values for the next one, but it doesn't seem to work.

Any clue of how to implement this or an example simple model implementation with such constraints would be of great help.

For the constraint |x(tk) - x(tk-1)| <= Cx for a state variable, create a new variable dx that is the defined as the derivative of the variable x and set a constraint on the ramp rate.

# State Variable
x = m.SV(value=0,ub=2)
dx = m.SV(value=0,lb=-0.5,ub=0.5)

# Process model
m.Equation(5*x.dt() == -x + 3*u)
m.Equation(dx == x.dt())

For the constraint 0.9 <= u(tk) / u(tk-1) <= 1.1 for the MV, set up a delay variable and constrain the value of u to be no more than a 10% change from the prior value.

u = m.MV(value=0.8, lb=0.1, ub=10) # avoid divide-by-zero
u.STATUS = 1  # allow optimizer to change

ud = m.Var(value=0.8)
m.delay(u,ud,1) # delay of 1
m.Equation(u<=1.1*ud)
m.Equation(u>=0.9*ud)

If the MV value u ever reaches zero, then the value won't be allowed to change further. It is generally safer to use something like DMAX or DMAXHI/DMAXLO to configure the rate of change of the MV. Below is a complete and minimal example that demonstrates the two constraint methods.

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

m = GEKKO()
m.time = np.linspace(0,10,21)

# Manipulated variable
u = m.MV(value=0.8, lb=0.1, ub=10) # avoid divide-by-zero
u.STATUS = 1  # allow optimizer to change
#u.DCOST = 0.1 # smooth out MV movement
#u.DMAX = 0.5   # limit MV rate of change

ud = m.Var(value=0.8)
m.delay(u,ud,1) # delay of 1
m.Equation(u<=1.1*ud)
m.Equation(u>=0.9*ud)

# State Variable
x = m.SV(value=0,ub=2)
dx = m.SV(value=0,lb=-0.5,ub=0.5)

# Process model
m.Equation(5*x.dt() == -x + 3*u)
m.Equation(dx == x.dt())

# Objective
m.Maximize(x)

m.options.IMODE = 6 # control
m.solve(disp=False)

plt.figure(figsize=(6,3.5))
plt.subplot(2,1,1)
plt.step(m.time,u.value,'b-',label='MV Optimized')
plt.legend(); plt.ylabel('Input')
plt.subplot(2,1,2)
plt.plot(m.time,x.value,'r--',label='SV Response')
plt.ylabel('State'); plt.xlabel('Time')
plt.legend(); plt.tight_layout()
plt.show()