I am trying to solve a BVP problem (Cosserat rod ODE) with gekko. The goal is to find the initial conditions nsol and msol (which correspond to the internal forces and moments of the rod) that minimize the cost function (the position of the final point of the rod), when integrating, the cosserat equations gives us P, R, nsol, msol, which correspond to the position, orientation, internal forces and moment in a section of the rod.
but I keep getting this error:
Exception: @error: Equation Definition Equation without an equality (=) or inequality (>,<) false STOPPING...
I am a beginner with gekko and although I have seen multiple threads with the same error, the source of the error seems to be different everytime. Could anyone please point me in the right direction ? Thank you very much
import numpy as np
import math
from scipy import integrate
import matplotlib.pyplot as plt
from gekko import GEKKO
E = 200e7
nu = 0.3
G = E/(2*(1+nu))
r = 0.01
rho = 8000
g = np.array([0, 0, 0])
ray = 1
A = np.pi*r**2
I = (np.pi*r**4)/4
J = 2*I
L = 1
Lfin = 1.5
Kse = np.diag([G*A, G*A, E*A])
Kbt = np.diag([E*I, E*I, G*J])
def antisym(y):
AS = np.array([[0, -y[2], y[1]], [y[2], 0, -y[0]], [-y[1], y[0], 0]])
return AS
m = GEKKO()
dl = 81
m.time = np.linspace(0, L, dl)
# Parameters
R = m.Array(m.Var, (3,3))
P = m.Array(m.Var, (3))
R[0,0].value = 1
R[1,1].value = 1
R[2,2].value = 1
R[0,1].value = 0
R[0,2].value = 0
R[1,0].value = 0
R[1,2].value = 0
R[2,0].value = 0
R[2,1].value = 0
P[0].value = 0
P[1].value = 0
P[2].value = 0
#R = m.Array(m.Var, (3,3),lb=0,ub=1, value = np.eye(3))
#P = m.Array(m.Var, (3), value = np.zeros(3))
v = m.Array(m.Var, (3))
u = m.Array(m.Var, (3))
# Variables
nsol = m.Array(m.Var, (3), value = 0)
msol = m.Array(m.Var, (3), value = 0)
test = np.zeros(dl)
test[-1] = 1.0
final = m.Param(value = test)
# Equations
m.Equation(v == np.dot(np.dot(np.diag((1/(G*A), 1/(G*A), 1/(E*A))), np.transpose(R)), nsol) + np.array([0,0,1]))
m.Equation(u == np.dot(np.dot(np.diag((1/(E*I), 1/(E*I), 1/(G*J))), np.transpose(R)), msol) + np.array([0,0,0]))
for i in range(2):
m.Equation(P[i].dt() == np.dot(R[i, :],v))
for i in range(2):
for j in range(2):
m.Equation(R[i, j].dt() == np.dot(R[i, :], antisym(u)[:, j]))
for i in range(2):
m.Equation(nsol[i].dt() == 0)
m.Equation(msol[0].dt() == -(P[1].dt()*nsol[2]-P[2].dt()*nsol[1]))
m.Equation(msol[1].dt() == -(P[2].dt()*nsol[0]-P[0].dt()*nsol[2]))
m.Equation(msol[2].dt() == -(P[0].dt()*nsol[1]-P[1].dt()*nsol[0]))
# Objective
m.Minimize(P[2]*final - Lfin)
m.options.IMODE = 6
m.solve()
One way to troubleshoot these types of errors is to inspect the gk0_model.apm model file in the run directory m.path. I modified the code to open the folder with m.open_folder() and the apm file:
Model
Parameters
p1
End Parameters
Variables
v1 = 1
v2 = 0
v3 = 0
v4 = 0
v5 = 1
v6 = 0
v7 = 0
v8 = 0
v9 = 1
v10 = 0
v11 = 0
v12 = 0
v13 = 0
v14 = 0
v15 = 0
v16 = 0
v17 = 0
v18 = 0
v19 = 0
v20 = 0
v21 = 0
v22 = 0
v23 = 0
v24 = 0
End Variables
Equations
False
False
$v10=((((v1)*(v13))+((v2)*(v14)))+((v3)*(v15)))
$v11=((((v4)*(v13))+((v5)*(v14)))+((v6)*(v15)))
$v1=((((v1)*(0))+((v2)*(v18)))+((v3)*((-v17))))
$v2=((((v1)*((-v18)))+((v2)*(0)))+((v3)*(v16)))
$v4=((((v4)*(0))+((v5)*(v18)))+((v6)*((-v17))))
$v5=((((v4)*((-v18)))+((v5)*(0)))+((v6)*(v16)))
$v19=0
$v20=0
$v22=(-((($v11)*(v21))-(($v12)*(v20))))
$v23=(-((($v12)*(v19))-(($v10)*(v21))))
$v24=(-((($v10)*(v20))-(($v11)*(v19))))
minimize (((v12)*(p1))-1.5)
End Equations
End Model
The first two equations are listed as False. This means that python evaluated the == is a comparative statement versus a symbolic expression. Gekko symbolic expressions are needed to compile the model into byte-code for automatic differentiation. In this case, the equations:
m.Equation(v == np.dot(np.dot(np.diag((1/(G*A), 1/(G*A), 1/(E*A))),\
np.transpose(R)), nsol) + np.array([0,0,1]))
m.Equation(u == np.dot(np.dot(np.diag((1/(E*I), 1/(E*I), 1/(G*J))),\
np.transpose(R)), msol) + np.array([0,0,0]))
are vectors and should be scalars.
# Equations
r1 = np.dot(np.dot(np.diag((1/(G*A), 1/(G*A), 1/(E*A))), \
np.transpose(R)), nsol) + np.array([0,0,1])
r2 = np.dot(np.dot(np.diag((1/(E*I), 1/(E*I), 1/(G*J))), \
np.transpose(R)), msol) + np.array([0,0,0])
for i in range(3):
m.Equation(v[i]==r1[i])
m.Equation(u[i]==r2[i])
This gives an unbounded solution error when attempting to solve. Adding a lower bound of -1000 and upper bound of 1000 for all variables gives a successful solution. If variables at at the bound, it may indicate that the problem is over-specified or unbounded without the artificial bounds.
import numpy as np
import math
from scipy import integrate
import matplotlib.pyplot as plt
from gekko import GEKKO
E = 200e7
nu = 0.3
G = E/(2*(1+nu))
r = 0.01
rho = 8000
g = np.array([0, 0, 0])
ray = 1
A = np.pi*r**2
I = (np.pi*r**4)/4
J = 2*I
L = 1
Lfin = 1.5
Kse = np.diag([G*A, G*A, E*A])
Kbt = np.diag([E*I, E*I, G*J])
def antisym(y):
AS = np.array([[0, -y[2], y[1]], [y[2], 0, -y[0]], [-y[1], y[0], 0]])
return AS
m = GEKKO()
dl = 81
m.time = np.linspace(0, L, dl)
# Parameters
R = m.Array(m.Var, (3,3), lb=-1000, ub=1000)
P = m.Array(m.Var, (3), lb=-1000, ub=1000)
R[0,0].value = 1
R[1,1].value = 1
R[2,2].value = 1
R[0,1].value = 0
R[0,2].value = 0
R[1,0].value = 0
R[1,2].value = 0
R[2,0].value = 0
R[2,1].value = 0
P[0].value = 0
P[1].value = 0
P[2].value = 0
#R = m.Array(m.Var, (3,3),lb=0,ub=1, value = np.eye(3))
#P = m.Array(m.Var, (3), value = np.zeros(3))
v = m.Array(m.Var, (3), lb=-1000, ub=1000)
u = m.Array(m.Var, (3), lb=-1000, ub=1000)
# Variables
nsol = m.Array(m.Var, (3), value = 0, lb=-1000, ub=1000)
msol = m.Array(m.Var, (3), value = 0, lb=-1000, ub=1000)
test = np.zeros(dl)
test[-1] = 1.0
final = m.Param(value = test)
# Equations
r1 = np.dot(np.dot(np.diag((1/(G*A), 1/(G*A), 1/(E*A))), \
np.transpose(R)), nsol) + np.array([0,0,1])
r2 = np.dot(np.dot(np.diag((1/(E*I), 1/(E*I), 1/(G*J))), \
np.transpose(R)), msol) + np.array([0,0,0])
for i in range(3):
m.Equation(v[i]==r1[i])
m.Equation(u[i]==r2[i])
for i in range(2):
m.Equation(P[i].dt() == np.dot(R[i, :],v))
for i in range(2):
for j in range(2):
m.Equation(R[i, j].dt() == np.dot(R[i, :], antisym(u)[:, j]))
for i in range(2):
m.Equation(nsol[i].dt() == 0)
m.Equation(msol[0].dt() == -(P[1].dt()*nsol[2]-P[2].dt()*nsol[1]))
m.Equation(msol[1].dt() == -(P[2].dt()*nsol[0]-P[0].dt()*nsol[2]))
m.Equation(msol[2].dt() == -(P[0].dt()*nsol[1]-P[1].dt()*nsol[0]))
# Objective
m.Minimize(P[2]*final - Lfin)
m.options.IMODE = 6
#m.open_folder()
m.solve()
Successful Solution Summary:
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 -1.2000000e+02 1.00e+00 1.24e-02 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -6.2000001e+02 4.70e-14 3.40e-01 -3.0 4.00e+04 - 6.60e-01 1.00e+00f 1
2 -1.1150000e+03 8.00e-14 6.43e-04 1.0 5.86e+04 - 1.00e+00 6.76e-01f 1
3 -1.1199121e+03 9.48e-14 3.86e-08 -1.1 3.93e+02 - 9.98e-01 1.00e+00f 1
4 -1.1199991e+03 7.96e-14 2.43e-10 -3.1 6.97e+00 - 9.98e-01 9.99e-01f 1
Reallocating memory for MA57: lfact (156431)
5 -1.1200000e+03 6.50e-14 2.43e-13 -9.0 7.03e-02 - 9.99e-01 9.99e-01f 1
Number of Iterations....: 5
(scaled) (unscaled)
Objective...............: -1.1200000091288521e+03 -1.1200000091288521e+03
Dual infeasibility......: 2.4264487412842937e-13 2.4264487412842937e-13
Constraint violation....: 6.4955110402786716e-14 6.4955110402786716e-14
Complementarity.........: 9.8229036600334927e-07 9.8229036600334927e-07
Overall NLP error.......: 9.8229036600334927e-07 9.8229036600334927e-07
Number of objective function evaluations = 6
Number of objective gradient evaluations = 6
Number of equality constraint evaluations = 6
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 6
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 5
Total CPU secs in IPOPT (w/o function evaluations) = 0.117
Total CPU secs in NLP function evaluations = 0.181
EXIT: Optimal Solution Found.
The solution was found.
The final value of the objective function is -1120.00000912885
---------------------------------------------------
Solver : IPOPT (v3.12)
Solution time : 0.334799999982351 sec
Objective : -1120.00000000000
Successful solution
---------------------------------------------------