Question: I have a system of differential equations that I'm trying to perform some optimal control on, using Gekko. In particular, I have a point-mass orbiting a planet and would simply like to raise its orbit using modelled thrusters as control inputs. In order to set the final radial position and velocity at the new raised orbit, I need to set the following boundary conditions: sqrt(x[-1]**2 + y[-1]**2) = r_final sqrt(vx[-1]**2 + vy[-1]**2) = v_final  Where (x,y) is the Cartesian position of the point-mass, and (vx,vy) is its velocity components, and the [-1] notation implies the last element in the trajectory. Furthermore, r_final is the desired altitude of the new orbit, and v_final is the desired speed of the new orbit. So far I have only been able to find functions in Gekko that constrain a single variable (e.g. fix(), fix_final() and fix_initial()), however I couldn't find any functions or examples that had more complex boundary conditions that included multiple variables. Looking at the above-mentioned example, I'm hoping to use Gekko to constrain the final radial position and speed at the new orbit similar to the following piece of code: m.fix(m.sqrt(x**2 + y**2), pos=len(m.time)-1,val=r_final) m.fix(m.sqrt(vx**2 + vy**2), pos=len(m.time)-1,val=v_final)  Is something like this possible? Thanks.

Answer: Only variables can be directly fixed in gekko. Below are three options. I recommend trying Option 2 first and then use Option 3 if the solver can't find a solution. Option 1: Define Equations With terminal complex expressions, use a final parameter that is non-zero only at the final time point with an equation definition. z = np.linspace(0,tf,n); z[-1]=1 final = m.Param(z) m.Equation(final*(x**2 + y**2 - r_final**2)==0) m.Equation(final*(m.sqrt(vx**2 + vy**2) - v_final)==0)  Don't use m.Equation(final*(x**2 + y**2)==r_final**2) or it will be infeasible with 0==r_final**2 before the final point. These equations create many 0==0 equations and can cause a problem for the gradient-based solvers thinking there are too few degrees of freedom. Option 2: Define Variables and fix_final() Another way to do this is to define two new variables and fix those at the final time point with fix_final(). r,v = m.Array(m.Var,2) m.Equations([r**2==x**2 + y**2, v**2==vx**2 + vy**2]) m.fix_final(r,r_final) m.fix_final(v,v_final)  Option 3: Minimize Soft Constraints If the solver isn't able to find a feasible solution, try setting a soft constraint instead of fix() or fix_final(). m.Minimize(final*100*(r-r_final)**2) m.Minimize(final*200*(v-v_final)**2)  Additional Info It is also possible to minimize final time for the maneuver by making tf a variable with m.Minimize(tf). The Bryson Benchmark problem (see Problem 2) has four examples of specifying terminal constraints.  from gekko import GEKKO import numpy as np import matplotlib.pyplot as plt #Build Model ############################ m = GEKKO() #define time space nt = 101 m.time = np.linspace(0,1.5,nt) #Parameters u = m.MV(value =0, lb = -1, ub = 1) u.STATUS = 1 u.DCOST = 0.0001 # slight penalty to discourage MV movement #Variables x1 = m.Var(value=0.5) x2 = m.Var(value = 0) myObj = m.Var() #Equations m.Equation(myObj.dt() == 0.5*x1**2) m.Equation(x1.dt() == u + x2) m.Equation(x2.dt() == -u) f = np.zeros(nt) f[-1] = 1 final = m.Param(value=f) # Four options for final constraints x1(tf)=0 and x2(tf)=0 option = 2 # best option = 3 (fix endpoints directly) if option == 1: # most likely to cause DOF issues because of many 0==0 equations m.Equation(final*x1 == 0) m.Equation(final*x2 == 0) elif option == 2: # inequality constraint approach is better but there are still # many inactive equations not at the endpoint m.Equation((final*x1)**2 <= 0) m.Equation((final*x2)**2 <= 0) elif option == 3: #requires GEKKO version >= 0.0.3a2 # fix the value just at the endpoint (best option) m.fix(x1,pos=nt-1,val=0) m.fix(x2,pos=nt-1,val=0) else: #penalty method ("soft constraint") that may influence the # optimal solution because there is just one combined objective # and it may interfere with minimizing myObj m.Obj(1000*(final*x1)**2) m.Obj(1000*(final*x2)**2) m.Obj(myObj*final) ######################################## #Set Global Options and Solve m.options.IMODE = 6 m.options.NODES = 3 m.options.MV_TYPE = 1 m.options.SOLVER = 1 # APOPT solver #Solve m.solve() # (remote=False) for local solution # Print objective value print('Objective (myObj): ' + str(myObj[-1])) ######################################## #Plot Results plt.figure() plt.plot(m.time, x1.value, 'y:', label = '$x_1$') plt.plot(m.time, x2.value, 'r--', label = '$x_2$') plt.plot(m.time, u.value, 'b-', label = 'u') plt.plot(m.time, myObj.value,'k-',label='Objective') plt.legend() plt.show()