Here is a similar problem with an inverted pendulum

There is also a related Stack Overflow question: how to use arrays in gekko optimizer for python Gekko can solve path planning problems that are subject to actuator constraints and equations of motion. One of the challenging mathematical features of this problem is how to smoothly pose the path constraints so that the robot does not pause at each intermediate path point. A potential method is to create a path cubic spline that helps to define the desired location for each time point. Here is one problem with matrices:
from gekko import GEKKO
import numpy as np
m = GEKKO(remote=False)
ni = 3; nj = 2; nk = 4
# solve AX=B
A = m.Array(m.Var,(ni,nj),lb=0)
X = m.Array(m.Var,(nj,nk),lb=0)
AX = np.dot(A,X)
B = m.Array(m.Var,(ni,nk),lb=0)
# equality constraints
m.Equations([AX[i,j]==B[i,j] for i in range(ni) \
for j in range(nk)])
m.Equation(5==m.sum([m.sum([A[i][j] for i in range(ni)]) \
for j in range(nj)]))
m.Equation(2==m.sum([m.sum([X[i][j] for i in range(nj)]) \
for j in range(nk)]))
# objective function
m.Minimize(m.sum([m.sum([B[i][j] for i in range(ni)]) \
for j in range(nk)]))
m.solve()
print(A)
print(X)
print(B)
Here is another test script that demonstrates dot products and the trace function with Numpy:
import numpy as np
from gekko import GEKKO
m = GEKKO(remote=False)
# Random 3x3
A = np.random.rand(3,3)
# Random 3x1
b = np.random.rand(3,1)
# Gekko array 3x3
p = m.Array(m.Param,(3,3))
# Gekko array 3x1
y = m.Array(m.Var,(3,1))
# Dot product of A p
x = np.dot(A,p) # or A@p
# Dot product of x y
w = x@y
# Dot product of p y
z = p@y # or np.dot(p,y)
# Trace (sum of diag) of p
t = np.trace(p)
# solve Ax = b
s = m.axb(A,b)
m.solve()