import numpy as np
from gekko import GEKKO

m = GEKKO()
m.time = np.linspace(0,20,41)

mass = 500
b = m.Param(value=50)
K = m.Param(value=0.8)

p = m.MV(value=0, lb=0, ub=100)
p.STATUS = 1  # allow optimizer to change
p.DCOST = 0.1 # smooth out gas pedal movement
p.DMAX = 20   # slow down change of gas pedal

v = m.CV(value=0)
v.STATUS = 1   # add the SP to the objective
m.options.CV_TYPE = 2 # sq error
v.SP = 40      # set point
v.TR_INIT = 1  # setpoint trajectory
v.TAU = 5      # time constant of setpoint trajectory

m.options.IMODE = 6 #control
m.Equation(mass*v.dt() == -v*b + K*b*p)
m.solve()

#%% Plot solution
import matplotlib.pyplot as plt
plt.figure()
plt.subplot(2,1,1)
plt.plot(m.time,p.value,'b-',lw=2)
plt.ylabel('gas')
plt.subplot(2,1,2)
plt.plot(m.time,v.value,'r--',lw=2)
plt.ylabel('velocity'); plt.xlabel('time')
plt.show()