title: How to use Gekko to solve for optimal control for a reusable reentry vehicle
body: <p>I am seeking to find optimal control (aoa and bank angle) to maximize cross range for a shuttle type reentry vehicle using Gekko.  Below is my code currently and I am getting a &quot;Solution not found&quot; with &quot;EXIT: Maximum Number of Iterations Exceeded&quot;.  The simulation assumes a point mass with a non-rotating earth frame. The EOMS are 6 coupled, non-linear ODEs. I have tried using different solvers, implementing/removing state and control constraints, increasing maximum number of iterations, etc.  I am not confident with my setup and implementation of the problem in Gekko and am hoping for some feedback on what I can try next. I have tried to follow the setup and layouts in  <a href="https://apmonitor.com/wiki/index.php/Main/GekkoPythonOptimization" rel="nofollow noreferrer">APMonitor's Example 11. Optimal Control with Integral Objective</a>, <a href="https://apmonitor.com/do/index.php/Main/InvertedPendulum" rel="nofollow noreferrer">Inverted Pendulum Optimal Control</a>, and <a href="https://apmonitor.com/wiki/index.php/Main/GekkoPythonOptimization" rel="nofollow noreferrer">Example 13. Optimal Control: Minimize Final Time</a>. Solutions I'm seeking are below.</p>
  <p>Any help is very much appreciated!</p>
  <pre><code>from gekko import GEKKO
  import numpy as np
  import matplotlib.pyplot as plt
  import math
  
  pi = math.pi
  
  ########################
  ######FRONT MATTER######
  ########################
  m = GEKKO() # initialize GEKKO
  nt = 2501 #simulation time is 2500 seconds
  tfin = 2500
  m.time = np.linspace(0,tfin,nt) #time array
  
  #==================#
  #PARAMS
  #==================#
  Re  = m.Param(value = 6371203.92)           # radius of the earth, m
  S   = m.Param(value = 249.9091776)          # vehicle surface area, m^2
  cl0 = m.Param(value = -0.2070)              # coeff lift param 1
  cl1 = m.Param(value = 1.6756)               # coeff lift param 2
  cd0 = m.Param(value = 0.0785)               # coeff drag param 1
  cd1 = m.Param(value = -0.3529)              # coeff drag param 2
  cd2 = m.Param(value = 2.0400)               # coeff drag param 3
  H   = m.Param(value = 7254.24)              # density scale height, m
  rho0= m.Param(value = 1.225570827014494)    # sea level atmospheric density, kg/m^3
  mu  = m.Param(value = 3.986031954093051e14) #earth gravitational param, m^3/s^2
  mass= m.Param(value = 92079.2525560557)     #vehicle mass, kg
  
  #===============================#
  #BOUNDARY CONDITIONS
  #===============================#
  t0     = 0
  alt0   = 79248
  rad0   = alt0+Re
  altf   = 24384
  radf   = altf+Re
  lon0   = 0
  lat0   = 0
  speed0 = +7802.88
  speedf = +762
  fpa0   = -1*pi/180
  fpaf   = -5*pi/180
  azi0   = +90*pi/180
  azif   = -90*pi/180
  
  #===============================#
  #LIMITS ON VARIABLES
  #===============================#
  tfMin = 0;              tfMax = 3000;
  radMin = Re;            radMax = rad0;
  lonMin = -pi;           lonMax = -lonMin;
  latMin = -70*pi/180;    latMax = -latMin;
  speedMin = 10;          speedMax = 45000;
  fpaMin = -80*pi/180;    fpaMax =  80*pi/180;
  aziMin = -180*pi/180;   aziMax =  180*pi/180;
  aoaMin = -90*pi/180;    aoaMax = -aoaMin;
  bankMin = -90*pi/180;   bankMax =   1*pi/180;
  
  #===============================#
  #VARIABLES
  #===============================#
  
  #state variables and bounds
  rad = m.Var(value=rad0, lb=radMin, ub=radMax)       # radius, m
  lon = m.Var(value=lon0, lb=lonMin, ub=lonMax)       # longitude, rad
  lat = m.Var(value=lat0, lb=latMin, ub=latMax)       # latitude, rad
  vel = m.Var(value=speed0, lb=speedMin, ub=speedMax) # velocity, m/sec
  fpa = m.Var(value=fpa0, lb=fpaMin, ub=fpaMax)       # flight path angle, rad
  azi = m.Var(value=azi0, lb=aziMin, ub=aziMax)       # azimuth angle, rad
  
  #control variables
  aoa     = m.MV(value=-20, lb=aoaMin, ub=aoaMax)    # angle of attack, rad
  aoa.STATUS = 1
  aoa.DCOST = 1e-2
  bank    = m.MV(value=0, lb=bankMin, ub=bankMax)    # bank angle, rad
  bank.STATUS = 1
  bank.DCOST = 1e-2
  
  #===============================#
  #INTERMEDIATE VARIABLES
  #===============================#
  altitude = m.Intermediate(rad - Re)
  CD       = m.Intermediate(cd0+cd1*aoa+cd2*aoa**2)
  rho      = m.Intermediate(rho0*m.exp(-altitude/H))
  CL       = m.Intermediate(cl0+cl1*aoa)
  q        = m.Intermediate(0.5*rho*vel**2)
  D        = m.Intermediate(q*S*CD/mass)
  L        = m.Intermediate(q*S*CL/mass)
  gravity  = m.Intermediate(mu/rad**2)
  
  #===============================#
  #EOMS
  #===============================#
  
  p = np.zeros(nt) # mark final time point
  p[-1] = 1.0
  final = m.Param(value=p)
  
  m.Equation(rad.dt() == vel*m.sin(fpa))
  m.Equation((rad*m.cos(lat))*lon.dt() == vel*m.cos(fpa)*m.sin(azi))
  m.Equation(rad*lat.dt() == vel*m.cos(fpa)*m.cos(azi))
  m.Equation(vel.dt() == -D-gravity*m.sin(fpa))
  m.Equation(vel*fpa.dt() == (L*m.cos(bank)-m.cos(fpa)*(gravity-vel**2/rad)))
  m.Equation(vel*azi.dt() == (L*m.sin(bank)/m.cos(fpa)+vel**2*m.cos(fpa)*m.sin(azi)*m.tan(lat)/rad))
  
  #===============================#
  #OPTIMIZATION SOLVER
  #===============================#
  m.Obj(-lat*final)
  
  
  # m.options.SOLVER = 3
  # m.options.IMODE = 6
  # m.solve(disp=True)
  
  m.options.MAX_ITER = 500
  m.options.IMODE = 6
  # m.options.NODES = 3
  # m.options.MV_TYPE = 1
  m.options.SOLVER = 3
  # m.open_folder()
  m.solve()
  
  </code></pre>
  <p>The solution should look something like this for altitude, velocity, AOA and bank angle: <a href="https://i.stack.imgur.com/ZpBnx.png" rel="nofollow noreferrer">Altitude</a> <a href="https://i.stack.imgur.com/4IuaZ.png" rel="nofollow noreferrer">Velocity</a> <a href="https://i.stack.imgur.com/9cc4R.png" rel="nofollow noreferrer">AngleofAttack</a> <a href="https://i.stack.imgur.com/XhMwf.png" rel="nofollow noreferrer">Bank</a></p>