Solving BVP problem with Gekko Error: @error: Equation Definition Equation without an equality (=) or inequality (>,<) false STOPPING
Solving BVP problem with Gekko Error: @error: Equation Definition Equation without an equality (=) or inequality (>,<) false STOPPING
我正在尝试用gekko解决BVP问题（Cosserat rod ODE）。
目标是找到使成本函数（杆的最终点的位置）最小化的初始条件 nsol 和 msol（对应于杆的内力和力矩），积分时，cosserat 方程给出 P , R, nsol, msol, 分别对应杆截面中的位置、方位、内力和力矩。
但我不断收到此错误：
Exception: @error: Equation Definition
Equation without an equality (=) or inequality (>,<)
false
STOPPING...
我是 gekko 的初学者，虽然我已经看到多个线程出现相同的错误，但错误的来源似乎每次都不同。
谁能指出我正确的方向？
非常感谢
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()
解决此类错误的一种方法是检查 运行 目录 m.path 中的 gk0_model.apm 模型文件。我修改了代码以使用 m.open_folder() 和 apm 文件打开文件夹：
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
前两个方程列为False。这意味着 python 评估了 == 是比较语句与符号表达式。需要 Gekko 符号表达式将模型编译成字节码以进行自动微分。在这种情况下，等式：
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]))
是向量，应该是标量。
# 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])
这会在尝试求解时出现无界求解错误。为所有变量添加 -1000 的下限和 1000 的上限可提供成功的解决方案。如果变量在边界上，则可能表明问题是过度指定的或没有人为边界的无界。
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()
成功解决方案总结：
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
---------------------------------------------------
我正在尝试用gekko解决BVP问题（Cosserat rod ODE）。 目标是找到使成本函数（杆的最终点的位置）最小化的初始条件 nsol 和 msol（对应于杆的内力和力矩），积分时，cosserat 方程给出 P , R, nsol, msol, 分别对应杆截面中的位置、方位、内力和力矩。
但我不断收到此错误：
Exception: @error: Equation Definition Equation without an equality (=) or inequality (>,<) false STOPPING...
我是 gekko 的初学者，虽然我已经看到多个线程出现相同的错误，但错误的来源似乎每次都不同。 谁能指出我正确的方向？ 非常感谢
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()
解决此类错误的一种方法是检查 运行 目录 m.path 中的 gk0_model.apm 模型文件。我修改了代码以使用 m.open_folder() 和 apm 文件打开文件夹：
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
前两个方程列为False。这意味着 python 评估了 == 是比较语句与符号表达式。需要 Gekko 符号表达式将模型编译成字节码以进行自动微分。在这种情况下，等式：
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]))
是向量，应该是标量。
# 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])
这会在尝试求解时出现无界求解错误。为所有变量添加 -1000 的下限和 1000 的上限可提供成功的解决方案。如果变量在边界上，则可能表明问题是过度指定的或没有人为边界的无界。
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()
成功解决方案总结：
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
---------------------------------------------------