GEKKO does not find optimal solution for a moon lander that does not go to zero height at zero speed

Question

I'm trying to solve an optimal control problem where a person lands on the Moon, but has a device that can propel her upwards, via a control parameter, alpha. The objective of the problem is to find the minimal time in which the moon lander can reach the surface of the moon, at speed zero (all motion along the vertical axis).

Now, I have implemented the code using gekko, with python, and it works just fine if the person starts, for instance, 40m above the surface of the Moon, and reaches the surface (final height = 0m) at speed zero. However, if I modify the code to have the person start, say, from 50m above the surface and get to a final height of 10m, gekko always converges to a point of local infeasibility. I have tried many different final heights, and it only seems to work when I set it to 0.

Is it a gekko problem or am I overlooking something in my code?

I followed the ideas shown here: https://apmonitor.com/do/index.php/Main/MinimizeFinalTime

Here's my code:

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO(remote=True)
m.solver_options = ['max_iter 500']
nt = 501
tm = np.linspace(0,1,nt)
m.time = tm

#variables

initialMass = m.Const(value=10.)
initialSpeed = m.Const(value=0.)
initialHeight = m.Const(value=40.)
finalHeight = m.Const(value=10.)

g = m.Const(value=1.625) #gravitational acceleration
k = m.Const(value=0.01)
power = m.Const(value=initialMass)


v = m.Var(value=initialSpeed) #position
h = m.Var(value=initialHeight) #height
mass = m.Var(value=initialMass, lb=0., ub=initialMass) #mass



#MV

alpha = m.MV(value=0.2, lb=0.2, ub=1.) #control variable
alpha.STATUS = 1

#FV

tf = m.FV(value=1., lb=0.1)
tf.STATUS = 1 # the value can be adjusted by the optimizer

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

#Equations

m.Equation(v.dt() == tf*(-g+2*alpha*g*power/mass))
m.Equation(h.dt() == tf*v)
m.Equation(mass.dt() == -tf*k*alpha)

m.Equation(h*final == finalHeight)
m.Equation(v*final == 0.)

m.Obj(tf)

m.options.IMODE = 6
m.solve(disp=True)
    

print('Solution')
print('Final time: '+str(tf.value[0]))

Answer 1

I managed to solve the problem. First, I added a lower bound for h, the height,

h = m.Var(value=initialHeight,lb=finalHeight)

and then I added an objective function to be minimized. Instead of

m.Equation(h*final == finalHeight)

now I have written

m.Minimize(final*(h+1)**2)

I followed the ideas presented here https://apmonitor.com/wiki/index.php/Apps/BrysonDenhamProblem