fix_final works for x_f=[0,0,0,0,0,0] but for absolutely no other final state - 'Solution Not Found'

'x_f' is defined in a separate file but here's the code for the GEKKO implementation:

m = GEKKO(remote = solverParams["remote"])
m.time = timeSeq
x = [m.Var(value = x_0[i], fixed_initial =  True) for i in range(len(x_0))]
u = [m.Var(value = u_0[i], fixed_initial = False) for i in range(len(u_0))]
w = m.Param(value = w)

for i in range(len(x)):
  m.fix_final(x[i], val = x_f[i])
  if stateBounds[i]["lower"] != "-Inf":
    x[i].lower = stateBounds[i]["lower"]
  if stateBounds[i]["upper"] != "+Inf":
    x[i].upper = stateBounds[i]["upper"]
for i in range(len(u)):
  if forceBounds[i]["lower"] != "-Inf":
    u[i].lower = forceBounds[i]["lower"]
  if forceBounds[i]["upper"] != "+Inf":
    u[i].upper = forceBounds[i]["upper"]
  
eqs = [x[i].dt() == np.matmul(A[i], x) + np.matmul(np.atleast_1d(B[i]), u)
        for i in range(len(x))]
eqs = m.Equations(eqs)

startTime = time.time()
m.Minimize(np.add(np.matmul(np.subtract(x, x_f), np.atleast_1d(np.matmul(np.atleast_1d(Q), np.subtract(x, x_f)))),
                  np.matmul(            u,       np.atleast_1d(np.matmul(np.atleast_1d(R),             u))))
          *w)
m.options.IMODE  = 6
m.options.SOLVER = 3
m.solve(disp = solverParams["disp"])
stopTime = time.time()

x_p = [x[i].value for i in range(len(x))]
u_p = [u[i].value for i in range(len(u))]
x_p = np.transpose(x_p)
u_p = np.transpose(u_p)

My test case is for:

x_0 = np.array([50.00, -25.00, 80.00, 0.00, 0.00, 0.00])

If I set x_f to be the exact zero vector, it gives me a solution but if I set it to, say,

x_f = np.array([ 0.04,   0.00,  0.00, 0.00, 0.00, 0.00])

it says that it hits a point of local infeasibility. I doubt that this is an actual scenario of infeasibility since this has happened for several cases of x_0 and x_f.

Here, the Q matrix is 0_{6x6} and R is 10*I_{3x3} and the system is linear time-invariant.

Fixing the final state also fixes the derivative of the final state to zero. This is a known issue with Gekko that has been reported elsewhere.

Instead of fixing the final value, try a terminal equation such as:

p = np.zeros_like(m.time); p[-1]=1
final = m.Param(p)
m.Equation(final*(x[i]-x_f[i])==0)

If this doesn't work, try a soft terminal constraint that is included in the objective function. This is typically the easiest for the solver and doesn't result in a potential infeasible solution.

m.Minimize(final*(x[i]-x_f[i])**2)

Here is an example of a double inverted pendulum that uses the soft terminal constraint to reach a final target.