Optimal control of constrained discrete Hamiltonian system. Using GEKKO?

Question

So in my master's thesis I am at the stage where I am to do a simple optimal control of my system. The dynamic system is of Hamiltonian form with constraints, and prior to looking at optimal control, the system has been solved stepwise using the RATTLE method, where it is solved at every time step by scipy's fsolve.

Now the idea is to add a control force, u, along with the partial derivative of the potential, U, and optimize on some objective function, to step through like this.

So I looked into some optimal control, and specifically GEKKO and have tried to understand how, or if it is possible to implement it. My first thought was to look at it as a dynamic control problem, but there isn't really any proper differential equation here? So I might just be slow or misunderstanding something, but is it possible to have a discrete evolution in the equations? If so, could someone explain a simple example?

Below is an example of my thought process:

import numpy as np
from gekko import GEKKO

m = GEKKO()

nt = 3
m.time = np.linspace(0,1,nt)

# Parameters
u = m.MV()
u.STATUS = 1
h = 1/(nt-1)

# Variables
y0 = m.Var(value=0)
yf = m.Var()

m.Equation(yf == y0 + h/2 *(2*y0 + u)) <---

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

# Objective Function
m.Obj(yf*final)

m.options.IMODE = 6
m.solve()

Is there some way to after the first step, update y0 to get the value of yf? I.e. the next step would be m.Equation(y2 == y1 + h/2 *(2*y1 + u)), and so on?

So that is the first question! Another thought was to not think of it as an evolution, but rather just make one set of constraint equations like this, and treat it as an optimization problem, not a dynamic control, but as the method is implicit, it doesn't directly spit out Y to be used as input in the next \Phi.

I feel like I am missing something simple here, does anyone have any thoughts or ideas?

A toy example that works, but seems unfeasible with any sort of real system:

from gekko import GEKKO

m = GEKKO()

# Parameters
u = m.MV(lb=-5, ub=5)
u.STATUS = 1
y0 = m.Param(value=0)

# Variables
y1 = m.Var()
y2 = m.Var()
y3 = m.Var()

m.Equation(y1 == y0 + u)
m.Equation(y2 == y1 + u)
m.Equation(y3 == y2 + u)

# Objective Function
m.Obj(y3)

m.options.IMODE = 3
m.solve()

Answer 1

Gekko has imode>=4 dynamic models such as DAEs (Differential and Algebraic equations), ARX, MPC (continuous or discrete), and other model forms. If they don't fit into one of these forms then there is the option to define the model variables and equations directly. This can be done efficiently with function definitions or arrays such as:

from gekko import GEKKO
m = GEKKO()
n = 4
u = m.Array(m.Var,n-1,lb=-5,ub=5)
y = m.Array(m.Var,n)
m.fix(y[0],0)
m.Equation(y[0]==0)
m.Equations([y[i+1] == y[i] + u[i] for i in range(n-1)])
m.Minimize((y[3]-2)**2)
m.options.IMODE = 3
m.solve()

There is no way to use imode>=4 and access each element of the time horizon in equations, but matrix and array operations make the model more compact.