Forming system of constraint equations using GEKKO optimization framework

I'm trying to solve a simple minimum time optimal control problem using double integrator dynamics of the form,

dx1/dt = x2
dx2/dt = u

with the GEKKO optimization framework as follows:

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

model = GEKKO(remote=False)

x1_initial = 0.0
x1_final = 10.0

x2_initial = 0.0
x2_final = 0.0

t_initial = 0.0
t_final = 25.0
num_timesteps = 1000
dt = (t_final - t_initial) / num_timesteps

x = model.Array(model.Var, (2, num_timesteps + 1))
u = model.Array(model.Var, num_timesteps + 1)
tf = model.Var()

for k in range(num_timesteps + 1):
    u[k].lower = -0.4
    u[k].upper = 0.4
    u[k].value = 0.0

for k in range(num_timesteps + 1):
    x[0, k].value = 5.0
    x[1, k].value = 0.0

tf.lower = t_initial
tf.upper = t_final
tf.value = t_final
dt = (tf - t_initial) / num_timesteps

def f(x, u, k):
    return np.array([x[1,k], u[k]])

for k in range(num_timesteps):
    model.Equations([x[:, k + 1] == x[:, k] + (dt/2.0)*(f(x, u, k + 1) + f(x, u, k))])
    # model.Equation(x[0, k + 1] == x[0, k] + (dt/2.0)*(x[1, k + 1] + x[1, k]))
    # model.Equation(x[1, k + 1] == x[1, k] + (dt/2.0)*(u[k + 1] + u[k]))

model.Equation(x[0, 0] == x1_initial)
model.Equation(x[0, num_timesteps] == x1_final)

model.Equation(x[1, 0] == x2_initial)
model.Equation(x[1, num_timesteps] == x2_final)

model.Minimize(tf)
model.options.solver = 3
model.solve()

# Plotting results
t = np.linspace(t_initial, tf.value, num_timesteps + 1)

u_optimal = []
for k in range(num_timesteps + 1):
    u_optimal.append(u[k].value)

x1_optimal = []
for k in range(num_timesteps + 1):
    x1_optimal.append(x[0, k].value)

x2_optimal = []
for k in range(num_timesteps + 1):
    x2_optimal.append(x[1, k].value)

plt.figure()
plt.plot(t, u_optimal)
plt.xlabel('time (s)')
plt.ylabel('u(t)')
plt.grid()

plt.figure()
plt.plot(t, x1_optimal)
plt.xlabel('time (s)')
plt.ylabel('x1(t)')
plt.grid()

plt.figure()
plt.plot(t, x2_optimal)
plt.xlabel('time (s)')
plt.ylabel('x2(t)')
plt.grid()

plt.show()

What I'm trying to do is to form a system of equality constraints using trapezoidal integration and then solve this system for the optimal control inputs using GEKKO. However, using the function definition,

def f(x, u, k):
    return np.array([x[1,k], u[k]])

in conjunction with the system of equality constraints,

for k in range(num_timesteps):
    model.Equations([x[:, k + 1] == x[:, k] + (dt/2.0)*(f(x, u, k + 1) + f(x, u, k))])

gives me the following error,

Exception: @error: Equation Definition
 Equation without an equality (=) or inequality (>,<)
 false
 STOPPING...

I've added two commented lines of code in the above code snippet that will allow the program to run correctly, but I'm hoping to avoid having to separate each equation out, since I'd like to extend this to problems that deal with more complicated system dynamics, and to also use more sophisticated collocation methods instead of the trapezoidal approach.

I know that GEKKO has some nice features for dynamic optimization, but I'm looking to try and implement various direct collocation methods myself to understand the theory a bit better.

BYU-PRISM / GEKKO

Forming system of constraint equations using GEKKO optimization framework #78