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
---------------------------------------------------
