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

# Initialize Model
m = GEKKO()

nt = 101
m.time = np.linspace(0,1,nt)

# Parameters
u = m.MV(value=-6)
u.STATUS = 1

# Variables
x = m.Var(value=0,ub=1/9)
v = m.Var(value=1)

J = m.Intermediate(0.5*m.integral(u**2))

p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)

# Equations
m.Equation(x.dt() == v)
m.Equation(v.dt() == u)

m.Equation(final*x==0)
m.Equation(final*(v+1)==0)

# ,or
# m.Minimize(final*1e5*x**2)
# m.Minimize(final*1e5*(v+1)**2)

# Objective Function
m.Minimize(final*J)

# Options
m.options.IMODE = 6
m.options.NODES = 2
m.options.SOLVER = 1

m.solve(disp=False)

print('Final Objective Function Value:', J.value[-1])

# Create a figure
plt.figure(figsize=(10,4))
plt.subplot(2,2,1)
plt.plot([0,1],[1/9,1/9],'r:',label=r'$x<\frac{1}{9}$')
plt.plot(m.time,x.value,'k-',lw=2,label=r'$x$')
plt.ylabel('Position')
plt.legend(loc='best')
plt.subplot(2,2,2)
plt.plot(m.time,v.value,'b--',lw=2,label=r'$v$')
plt.ylabel('Velocity')
plt.legend(loc='best')
plt.subplot(2,2,3)
plt.plot(m.time,u.value,'r--',lw=2,label=r'$u$')
plt.ylabel('Thrust')
plt.legend(loc='best')
plt.xlabel('Time')
plt.subplot(2,2,4)
plt.plot(m.time,J.value,'g-',lw=2,label=r'$\frac{1}{2} \int u^2$')
plt.text(0.5,3.0,'Final Value = '+str(np.round(J.value[-1],2)))
plt.ylabel('Objective')
plt.legend(loc='best')
plt.xlabel('Time')
plt.show()