[Python] I want to simulate linear MPC control of a FOPTD model using GEKKO

I'm simulating a first-order plus time delay (FOPTD) process and controlling it with a GEKKO MPC model using identical parameters for both the plant and model (i.e. no model mismatch). However, when a nonzero time delay is introduced, the simulation becomes unstable.

• The process is modeled as a FOPTD system with gain K=3, time constant tau=7, and delay theta=5.
• The MPC model in GEKKO uses the same dynamics.
• I'm using Euler integration to simulate the process.
• GEKKO's m.delay() is used to model the input delay in the MPC.

Full Minimal Reproducible Example (MWE)

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

# Simulated process (nonlinear plant)
def process(t, t_list, u_list, y):
u0 = u_list[0]
y0 = 5.0
K = 3.0
tau = 7.0
theta = 5.0 # seconds

delay_time = t - theta
u = u0 if delay_time < 0 else np.interp(delay_time, t_list, u_list)

y = y + dt * ((-(y - y0) + K * (u - u0)) / tau)
return y

# GEKKO MPC model (same parameters)
K_model = 3.0
tau_model = 7.0
theta_model = 5.0 # same as plant delay

m = GEKKO(remote=False)
tmax_gekko = 10
dt_gekko = 1
m.time = np.linspace(0, tmax_gekko, int(tmax_gekko/dt_gekko) + 1)

m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 2 # trajectory tracking

u_gekko = m.MV(2)
u_gekko.STATUS = 1

ud = m.Var() # delayed input

y_gekko = m.CV(5)
y_gekko.STATUS = 1
y_gekko.FSTATUS = 1

# Apply delay
m.delay(u_gekko, ud, int(theta_model / dt_gekko) + 1)

# FOPTD model equation
m.Equation(y_gekko.dt() == (-(y_gekko - 5.0) + K_model*(ud - 2.0)) / tau_model)

# Simulation
tmax = 50
dt = 1
nt = int(tmax/dt) + 1
tspan = np.linspace(0, tmax, nt)

y_list = [5.0] + [0]*(nt-1)
u_list = [0]*nt

for i in range(nt - 1):
t = tspan
y_gekko.MEAS = y_list
y_gekko.SP = 5.0

m.solve(disp=False)
u = u_gekko.NEWVAL
u_list = u

y = process(t, tspan[:i+1], u_list[:i+1], y_list)
y_list[i+1] = y

print(f"Progress: {100*i/nt:.1f}%", end='\r')

u_list[-1] = u_list[-2]

# Plot results
fig, axs = plt.subplots(2, 1, sharex=True)
axs[0].plot(tspan, y_list, label='Process Output y')
axs[0].legend()
axs[1].plot(tspan, u_list, 'r', label='Control Input u')
axs[1].legend()
plt.tight_layout()
plt.show()


Why does the controller cause instability when a nonzero time delay is introduced, even though the process and model are perfectly matched?

• Is m.delay() not working as expected in this case?
• Do I need to modify the MV or CV tuning or add constraints?
• Does GEKKO handle integer time delays reliably in this simulation setup?

Continue reading...