# prob4
# Using GEKKO to solve Optimal Control Problems in MISER Manual (Problem 4)
# 
# we are trying to introduce the use of GEKKO to native MISER user
# this file is prepared by Dr. Joseph Lee
# 
# Dr. H. W. J. Lee
# Department of Applied Mathematics
# The Hong Kong Polytechnic University

# if gekko is not installed, uncomment (delete #) in the following line
#! pip install gekko

# only need to do it once, and put it (#) back as comment once gekko is 
# installed, e.g for anaconda python, you only need to do it once to 
# install gekko
#
# however, if you are using Google CoLab to run python, 
# then, you need to install it every time, if so, keep uncomment if you are
# using Google CoLab

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

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

z1 = m.FV(value=1,lb=1.0,ub=15.0)
z1.STATUS = 1
z2 = m.FV(value=0.0,lb=0.5,ub=10.0)
z2.STATUS = 1


x1 = m.Var(value=0)
x2 = m.Var(value=1)
x3 = m.Var(value=0)
x4 = m.Var(value=0)
myObj = m.Var()


u = m.MV(value=0,lb=-10,ub=10)
u.STATUS = 1
m.Equation(myObj.dt() == -z1)
m.Equation(x1.dt()==  x2              )
m.Equation(x2.dt()== -z1*x1/((x3+z2)**2) )
m.Equation(x3.dt()==  u               )
m.Equation(x4.dt()==  x3+z2           )

m.Equation(x3+z2-0.5 >=0)

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

m.Equation(x1*final == 0)
m.Equation((x4-1)*final == 0)

m.Obj(myObj)
m.options.IMODE = 6

m.solve(disp=False)

print('z1=',z1.value[-1],'   z2=',z2.value[-1])

z1= 12.770310983    z2= 0.50000924217

plt.figure(figsize=(15,8))
plt.plot(m.time, x1.value, 'y:', label = '$x_1$')
plt.plot(m.time, x2.value, 'r--', label = '$x_2$')
x3z=np.zeros(len(m.time))
constraint=np.zeros(len(m.time))
for number in range(len(m.time)):
    x3z[number]=x3.value[number]+z2.value[number]
    constraint[number]=0.5
plt.plot(m.time, x3z, 'g--', label = '$x_3+z2$')
plt.plot(m.time, x4.value, 'b--', label = '$x_4$')
plt.plot(m.time, constraint, 'm--', label = '$0.5$')

plt.plot(1.0,0.0,'ro') 
plt.plot(1.0,1.0,'ro') 

plt.legend()

plt.show()

plt.figure(figsize=(15,8))
plt.plot(m.time, u.value, 'b-', label = 'u')
plt.legend()

plt.show()