Gekko optimal control. I can't keep one value within its bounds

I am trying to simulate 3D plane flight. I have an issue with gamma value (Flight-path angle). It gets out of its bounds and then, the simulation stops. The gamma value is being calculated by this equation:

I turned it into this: m.Equation(gamma.dt()==tf*((L*m.cos(Mu)-mass*g*m.cos(gamma))/mass*V))

The target of the simulation is for the plane to reach certain X an Y values(m.Minimize(w*final*(x-pathx)**2) and m.Minimize(w*final*(pathy-y)**2)), while minimizing fuel consumed m.Maximize(0.2*mass*tf*final).

The solver controls gamma value, by controlling lift coefficient Cl, which affects the lift value L, which, in turn, affects gamma value. The equation that calculates lift L value looks like this: m.Equation(L==0.5*Ro*(V**2)*(Cl)*S). But in the end the solver can not control gamma value till the plane gets to its destination.

What could be messing with it?

My code:

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
import math
#Gekko model
m = GEKKO(remote=False)

#Time points
nt = 11
tm =  np.linspace(0,100,nt)
m.time = tm

# Variables
Ro = m.Var(value=1.1)#air density
g = m.Const(value=9.80665)
pressure = m.Var(value=101325)#
T = m.Var(value=281,lb=100)#temperature
T0 = m.Const(value=288)#temperature at see level
S = m.Const(value=122.6)
Cd = m.Var(value=0.025)#drag coef 0.06 works
#Cl = m.Const(value=0.3)#lift couef
FuelFlow = m.Var()
D = m.Var(value=10000,lb=0)#drag
Thrmax = m.Const(value=200000)#maximum throttle
Thr = m.Var()
V = m.Var(value=200,lb=45,ub=240)#velocity
gamma = m.Var(value=0,lb=-0.6,ub=1.2)# Flight-path angle
#gammaa = gamma.value
Xi = m.Var(value=0, lb=-2, ub=2.0)# Heading angle
x = m.Var(value=0,lb=-1000,ub=1015000)#x position
y = m.Var(value=0,lb=-1000,ub=1011000)#y position
h = m.Var(value=1000,lb=-20000,ub=50000)# height
targetH = m.Param(value=10000) #target flight altitude
mass = m.Var(value=60000,lb=10000)
pathx = m.Const(value=1000000) #intended distance in x direction
pathy = m.Const(value=1000000) #intended distance in y direction
L = m.Var(value=250000)#lift

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

m.options.MAX_ITER=10000 # iteration number

#Fixed Variable
tf = m.FV(value=1,lb=0.0001,ub=100.0)#
tf.STATUS = 1

# Controlled parameters
Tcontr = m.MV(value=0.6,lb=0.0,ub=1)# solver controls throttle pedal position
Tcontr.STATUS = 1
Tcontr.DCOST = 0

Mu = m.MV(value=0,lb=-1,ub=1)# solver controls bank angle 
Mu.STATUS = 1
Mu.DCOST = 1e-3

Cl = m.MV(value=0.1,lb=-0.3,ub=0.9)# solver controls lift couef 
Cl.STATUS = 1
Cl.DCOST = 1e-3

# Equations
m.Equation(Thr==Tcontr*Thrmax)
m.Equation(FuelFlow==0.75882*(1+(V/2938.5)))
m.Equation(D==0.5*Ro*(V**2)*Cd*S)
m.Equation(mass.dt()==tf*(-Thr*(FuelFlow/60000)))#
m.Equation(V.dt()==tf*((Thr-D-mass*g*m.cos(gamma))/mass)) # 
m.Equation(x.dt()==tf*(V*(m.cos(gamma))*(m.cos(Xi))))#
m.Equation(x*final<=pathx)
#pressure and density part
m.Equation(T==T0-(0.0065*h))
m.Equation(pressure==101325*(1-(0.0065*h)/T0)**((g*0.0289652)/(8.31446*0.0065)))#
m.Equation(Ro*(8.31446*T)==(pressure*0.0289652))
#2D addition part
m.Equation(L==0.5*Ro*(V**2)*(Cl)*S)# Problem here or no problem, idk
m.Equation(mass*m.cos(gamma)*V*Xi.dt()==tf*((L*m.sin(Mu)))) # 
m.Equation(y.dt()==tf*(V*(m.cos(gamma))*(m.sin(Xi))))#
#m.Equation((y-pathy)*final==0) 
m.Equation(y*final<=pathy)
#3D addition part
m.Equation(h.dt()==tf*(V*(m.sin(gamma))))#
m.Equation(h*final<=targetH)
m.Equation(gamma.dt()==tf*((L*m.cos(Mu)-mass*g*m.cos(gamma))/mass*V))#
#Cd equation
m.Equation(Cd==((Cl)**2)/10)


# Objective Function
w = 1e4
m.Minimize(w*final*(x-pathx)**2) #1D part (x)
m.Minimize(w*final*(pathy-y)**2) #2D part (y)

m.Maximize(0.2*mass*tf*final) #objective function
m.options.IMODE = 6
m.options.NODES = 2 # it was 3 before
m.options.MV_TYPE = 1
m.options.SOLVER = 3
#m.open_folder() # to search for infeasibilities
m.solve()
    
tm = tm * tf.value[0]
fig, axs = plt.subplots(11)
fig.suptitle('Results')
axs[0].plot(tm,Tcontr,'r-',LineWidth=2,label=r'$Tcontr$')
axs[0].legend(loc='best')
axs[1].plot(tm,V.value,'b-',LineWidth=2,label=r'$V$')
axs[1].legend(loc='best')
axs[2].plot(tm,x.value,'r--',LineWidth=2,label=r'$x$')
axs[2].legend(loc='best')
axs[3].plot(tm,D.value,'g-',LineWidth=2,label=r'$D$')
axs[3].legend(loc='best')
axs[4].plot(tm,L.value,'p-',LineWidth=2,label=r'$L$')
axs[4].legend(loc='best')
axs[5].plot(tm,y.value,'p-',LineWidth=2,label=r'$y$')
axs[5].legend(loc='best')
axs[6].plot(tm,Ro.value,'p-',LineWidth=2,label=r'$Ro$')
axs[6].legend(loc='best')
axs[7].plot(tm,h.value,'p-',LineWidth=2,label=r'$h$')
axs[7].legend(loc='best')
axs[8].plot(tm,gamma.value,'p-',LineWidth=2,label=r'$gamma$')
axs[8].legend(loc='best')
axs[9].plot(tm,Cl.value,'p-',LineWidth=2,label=r'$Cl$')
axs[9].legend(loc='best')
axs[10].plot(tm,Cd.value,'p-',LineWidth=2,label=r'$Cd$')
axs[10].legend(loc='best')    
plt.xlabel('Time')
#plt.ylabel('Value')
plt.show()

Important update

I want the solver to be able to keep gamma value within its boundaries till the necessary x and y values are achieved. The solver has issues with that.

Exhibit A.

In this case the simulation stopped because it could not keep gamma from exceeding its lower bound. gamma usually gets bigger when L (lift) * cos(mu) is bigger than mass*g*cos(gamma), and its gets lower when the situation is the opposite: mass*g*cos(gamma) is bigger than the L (lift) * cos(mu). Then the question becomes: why did suddenly mass*g*cos(gamma) became so much bigger than the L (lift) * cos(mu)? The g was constant. The mass did not change that much during the last moments of the simulation. The L (lift) did not become particularly small. cos(mu) is usually equal to 1 during this part of the simulation.

Exhibit B. In this case the simulation stopped, once again, because it could not keep gamma value within its upper bound. It is visible , that Cl value during the last moments of the simulation was rising, it was necessary to keep gamma from exceeding its lower bound, but after that, the gamma value spiked and, for some reason, the solver did not lower Cl value, which forced the gamma value to exceed its upper bound, which forced the simulation to stop before the target goals are achieved.

In this case the issue is: Why does doesn't the solver lower the Cl value to stop gamma from exceeding its upper bound?

from Gekko optimal control. I can't keep one value within its bounds