GEKKO Dynamic Optimisation Infeasible below certain timestep length

I want to do dynamic simulation of Electrocoagulation Process controlled with PI Controller in order to test whether my model has worked well and can be used for further Optimization Problem that I want to solve.

The MV is Electrical Current, which I write down as ordinary Variable m.Var. The CV is ADC Reading of Turbidity Sensor, which I also write down as an ordinary Variable m.Var.

Since the DOF is negative, I couldn't use IMODE = 4, so I use IMODE = 6 instead but without Objective Function. Also, I use APOPT here.

I tried at least 2 scenarios to test my model:

1. Current with no upper bound, total timestep nt = 501, with various final time t_final (from 2000 - 50000 seconds). The goal is to see when Turbidity attain Steady-State.
2. Current with upper bound of 5.2 amperes, total timestep nt = 501, with various final time (from 2000 - 300000 seconds). The goal is the same, but the upper bound represent the operational limit of the plant.

Scenario 1 proceed seamlessly. However, during Scenario 2, with final time 2000 s, the solver said it is infeasible. I initially thought it was due to the Turbidity lower bound being violated, thus I removed it. Yet, what I found afterwards are as follows:

1. Infeasibility persisted after Turbidity lb removed. (It doesn't cause infeasibility)
2. Problem Solved only when time interval is sparse (e.g. nt = 41 with t_final = 2000 s, or nt = 501 but with t_final = 300000 s)

When I tried to refine the time interval by increasing nt or decreasing t_final, the solver always fail to solve the problem.

How could this happen? What are the explanation behind this behavior? and What should I do to resolve the issue?

Addition: The sample code I am working now is this.

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

# Model GEKKO
m = GEKKO(remote=False)

# Constants and Parameters
## Time
# (I changed nt and tf repeatedly to see their dynamic characteristics, error, etc.)
nt = 11   # Data Points
tf = 5000   # [sec] Final Time
m.time = np.linspace(0,tf,nt)   # Time Array
t = m.Param(value=m.time)       # Time Parameter

## Electrocoagulations
Turb_ref = m.Param(value=100)  # [NTU] Turbidity Setpoint
Turb0 = 800                    # [NTU] Initial Value of Turbidity
ADC_ref = (Turb_ref-1323)/(-1.8457) # [ADC] ADC Value of Turbidity Sensor when Turbidity = 100 NTU
ADC0 = (Turb0-1323)/(-1.8457)  # [ADC] Initial Value of the ADC Turbidity Sensor Readings
CurrentMin = 1.8  # [Ampere] Minimum Electric Current
Current0 = 2.5 # [Ampere] Initial Current

# FOPDT Parameter
Kp = 909.8     # Process Gain
taup = m.Param(value=844.7)  # Process Time Constant
thetap = 395.7 # Process Dead Time

#################################################################
# Variables
## Elektrocoagulation Dynamic Process
# (I tested whether upper bound of Current affecting its dynamic behaviour or sim. error, etc.)
Current = m.Var(value=Current0,lb=1.8,ub=5.2)  # [Ampere] MV(Input): Electric Current
#Current = m.Var(value=Current0,lb=1.8)  # [Ampere] MV(Input): Electric Current - No Upper Bound
ADC = m.Var(value=ADC0,lb=0)  # [ADC] CV(Output): Turbidity Sensor Reading in ADC
Turb = m.Var(value=Turb0)   # [NTU] Turbidity for us to observe

# Controller
Kc = m.FV(value=8.8183e-2)    # Controller Gain
tauc = m.FV(value=18980)  # Controller Integral Time
Kc.STATUS = 0           # (I want to make the value of Controller Parameters constant the whole time)
tauc.STATUS = 0
ie = m.Var()    # Integral Error
e = m.Var()     # Error

#################################################################
# Equations

# Intermediate for FOPDT Dead Time Delay
alive = m.Intermediate(1 - m.exp(-(t-thetap)/taup))
# Process Dynamics
m.Equations([
    ADC == ADC0 + Kp * (Current - CurrentMin) * m.if3(t-thetap,0,alive),  # Electrocoagulation Process
    Turb == -1.8457*ADC + 1323  # Convert ADC to Turbidity for User to Observe [NTU]
    ])
# Feedback Control
m.Equations([
    e == ADC_ref - ADC,   # Error
    ie.dt() == e,      # Integration Error
    Current == Current0 + Kc*(e + ie/tauc),  # Control Signal Output
])

#################################################################

# Objective
#m.Obj(Arus)

# Solver Options
m.options.SOLVER = 1
m.options.IMODE = 6
#m.options.NODES = 4
#m.options.SENSITIVITY = 1
#m.open_folder() # open folder if remote=False to see infeasibilities.txt

# Solver
m.solve(GUI=True)