3

我知道这听起来很简单,但我正在为一个应该可以工作的高中课程编写一些代码。基本上,我正在解决一个非常受限的定义金字塔的方程组。问题是求解器说我的一个方程没有一个等式,尽管看起来它们都有。

我已经尝试了几乎所有的东西。我重新输入了方程式并仔细检查了所有语法。我已经有 2 天了同样的错误,我非常卡住

#sorry that the equations are messy but, they all have only one equality.

from gekko import GEKKO
V = 20
S = 30
t1 = 0.1
t2 = 0.2
t3 = 0.3
t4 = 0.05
t5 = 0.6
t6 = 0.2
t7 = 0.1


m = GEKKO()
x1= m.Var(value = 0)
x2= m.Var(value = 0)
x3= m.Var(value = 0)
y1= m.Var(value = 0)
y2= m.Var(value = 0)
y3= m.Var(value = 0)
z1= m.Var(value = 0)
z2= m.Var(value = 0)
z3= m.Var(value = 0)
t8= m.Var(value = 0)
t9= m.Var(value = 0)
t10= m.Var(value = 0)
t11= m.Var(value = 0)
t12= m.Var(value = 0)
m.Equations([
            V == z3*(x1*y2 - x2*y1) + y3*(x2*z1 - x1*z3) + x3*(y1*z3 - y2*z1),\
            S == 0.5 * m.sqrt(((x1*y2-x2*y1)**2)+((x2*z1-x1*z2)**2)+((x1*y2-x2*y1)**2)) + 0.5 * m.sqrt(((x1*y3-x3*y1)**2)+((x3*z1-x1*z3)**2)+((y1*z3-y3*z1)**2))+ 0.5 * m.sqrt(((x1*y3-x3*y1)**2)+((x3*z1-x1*z3)**2)+((y1*z3-y3*z1)**2)),\
            t1 == [m.acos(x1*x2+y1*y2+z1*z2)/((m.sqrt((x1**2) +(y1**2) + (z1**2)))*(m.sqrt((x2**2) +(y2**2) + (z2**2))))],\
            t2 == [m.acos(x1*x3+y1*y3+z1*z3)/((m.sqrt((x1**2) +(y1**2) + (z1**2)))*(m.sqrt((x3**2) +(y3**2) + (z3**2))))],\
            t3 == [m.acos(x2*x3+y2*y3+z2*z3)/((m.sqrt((x2**2) +(y2**2) + (z2**2)))*(m.sqrt((x3**2) +(y3**2) + (z3**2))))],\
            t4 == [m.acos((x1*(x1-x3)) +(y1*(y1-y3)) +(z1*(z1-z3))/((m.sqrt((x1**2) + (y1**2)+(z1**2))) * (m.sqrt((x1-x3)**2  + (y1-y3)**2 + (z1-z3)**2))))],\
            t5 == [m.acos((x1*(x1-x2))+ (y1*(y1-y2))+(z1*(z1-z2))/((m.sqrt((x1**2) + (y1**2) + (z1**2)))*(m.sqrt((x1-x2)**2 + (y1-y2)**2 +(z1-z2)**2))))],\
            t6 == [m.acos((x2*(x1-x2))+ (y2*(y1-y2))+(z2*(z1-z2))/((m.sqrt((x2**2) +(y2**2) + (z2**2)))*(m.sqrt((x1-x2)**2 +(y1-y2)**2 +(z1-z2)**2))))],\
            t7 == [m.acos((x2*(x2-x3))+ (y2*(y2-y3))+(z2*(z2-z3))/((m.sqrt((x2**2) +(y2**2) + (z2**2)))*(m.sqrt((x2-x3)**2 +(y2-y3)**2 +(z2-z3)**2))))],\
            t8 == [m.acos((x3*(x2-x3))+ (y3*(y2-y3))+(z3*(z2-z3))/((m.sqrt((x3**2) +(y3**2) + (z3**2)))*(m.sqrt((x2-x3)**2 +(y2-y3)**2 +(z2-z3)**2))))],\
            t9 == [m.acos((x1*(x1-x3))+ (y1*(y1-y3))+(z1*(z1-z3))/((m.sqrt((x1**2) +(y1**2) + (z1**2)))*(m.sqrt((x1-x3)**2 +(y1-y3)**2 +(z1-z3)**2))))],\
            t10 == [m.acos(((x1-x3) * (x2-x3) + (y1-y3) * (y2-y3) + (z1-z3) * (z2-z3)) /(m.sqrt((x1-x3)**2 + (y1-y3)**2 + (z1-z3)**2) * m.sqrt((x2-x3)**2 +(y2-y3)**2 + (z2-z3)**2)))],\
            t11 == [m.acos(((x1-x3)* (x1-x2) + (y1-y3)*(y1-y2) +(z1-z3)*(z1-z2))/(m.sqrt((x1-x3)**2 + (y1-y3)**2 + (z1-z3)**2)* m.sqrt((x1-x2)**2 +(y1-y2)**2 +(z1-z2)**2)))],\
            t12 == [m.acos(((x1-x2) *(x2-x3)+(y1-y2)*(y2-y3) +(z1-z2)*(z2-z3))/(m.sqrt((x1-x2)**2 +(y1-y2)**2 + (z1-z2)**2)* m.sqrt((x2-x3)**2 +(y2-y3)**2 +(z2-z3)**2)))] ])
m.solve(disp = True)
print("vector a:",x1.value,y1.value,z1.value, " vector b:",x2.value,y2.value,z2.value," vector c:",x3.value,y4.value,z3.value,)

错误是:异常:@error:方程式定义方程式没有等式(=)或不等式(>,<)false

但我不知道哪个方程不满足要求。

4

1 回答 1

0

问题是方程中的额外括号:

t1 == [m.acos(....)]

您可以通过删除等式右侧定义列表的括号来解决此问题。

t1 == m.acos(....)

下面是一个修改后的脚本,它仍然返回一个不可行的解决方案,但有一些额外的东西可以帮助您诊断问题。

  • 命名变量,以便 infeasibilities.txt 文件更具可读性
  • 用 \ 包裹方程,使其更具可读性
  • 检查所有方程,例如第二个方程,第 3 部分。应该是+ 0.5 * m.sqrt(((x2*y3-x3*y2)**2)吗?它目前看起来像是该等式第 2 部分的副本。
  • lb用(下限)和ub(上限)对变量设置现实的界限。例如,如果它是一个距离,那么下限可能是 0.01 或其他一些小数。
  • 包括更好的初始猜测值。从零开始可能意味着方程最初是在无穷大或 NaN 处计算的。
from gekko import GEKKO
V = 20
S = 30
t1 = 0.1
t2 = 0.2
t3 = 0.3
t4 = 0.05
t5 = 0.6
t6 = 0.2
t7 = 0.1

m = GEKKO(remote=False)
x1= m.Var(value = 0.11, lb=0.01, ub=2.0, name='x1')
x2= m.Var(value = 0.12, lb=0.01, ub=2.0, name='x2')
x3= m.Var(value = 0.13, lb=0.01, ub=2.0, name='x3')
y1= m.Var(value = 0.14, lb=0.01, ub=2.0, name='y1')
y2= m.Var(value = 0.15, lb=0.01, ub=2.0, name='y2')
y3= m.Var(value = 0.16, lb=0.01, ub=2.0, name='y3')
z1= m.Var(value = 0.17, lb=0.01, ub=2.0, name='z1')
z2= m.Var(value = 0.18, lb=0.01, ub=2.0, name='z2')
z3= m.Var(value = 0.19, lb=0.01, ub=2.0, name='z3')
t8= m.Var(value = 0.20, lb=0.01, ub=2.0, name='t8')
t9= m.Var(value = 0.21, lb=0.01, ub=2.0, name='t9')
t10= m.Var(value = 0.22, lb=0.01, ub=2.0, name='t10')
t11= m.Var(value = 0.23, lb=0.01, ub=2.0, name='t11')
t12= m.Var(value = 0.24, lb=0.01, ub=2.0, name='t12')
m.Equations([
            V == z3*(x1*y2 - x2*y1) + y3*(x2*z1 - x1*z3) + x3*(y1*z3 - y2*z1),\
            S == 0.5 * m.sqrt(((x1*y2-x2*y1)**2)+((x2*z1-x1*z2)**2)+((x1*y2-x2*y1)**2)) \
               + 0.5 * m.sqrt(((x1*y3-x3*y1)**2)+((x3*z1-x1*z3)**2)+((y1*z3-y3*z1)**2)) \
               + 0.5 * m.sqrt(((x1*y3-x3*y1)**2)+((x3*z1-x1*z3)**2)+((y1*z3-y3*z1)**2)),\
            t1 == m.acos(x1*x2+y1*y2+z1*z2)/((m.sqrt((x1**2) +(y1**2) \
                  + (z1**2)))*(m.sqrt((x2**2) +(y2**2) + (z2**2)))),\
            t2 == m.acos(x1*x3+y1*y3+z1*z3)/((m.sqrt((x1**2) +(y1**2) \
                  + (z1**2)))*(m.sqrt((x3**2) +(y3**2) + (z3**2)))),\
            t3 == m.acos(x2*x3+y2*y3+z2*z3)/((m.sqrt((x2**2) +(y2**2) \
                  + (z2**2)))*(m.sqrt((x3**2) +(y3**2) + (z3**2)))),\
            t4 == m.acos((x1*(x1-x3)) +(y1*(y1-y3)) +(z1*(z1-z3))/  \
                  ((m.sqrt((x1**2) + (y1**2)+(z1**2))) * (m.sqrt((x1-x3)**2  \
                  + (y1-y3)**2 + (z1-z3)**2)))),\
            t5 == m.acos((x1*(x1-x2))+ (y1*(y1-y2))+(z1*(z1-z2))/ \
                  ((m.sqrt((x1**2) + (y1**2) + (z1**2)))*(m.sqrt((x1-x2)**2 \
                  + (y1-y2)**2 +(z1-z2)**2)))),\
            t6 == m.acos((x2*(x1-x2))+ (y2*(y1-y2))+(z2*(z1-z2))/ \
                  ((m.sqrt((x2**2) +(y2**2) + (z2**2)))*(m.sqrt((x1-x2)**2 \
                  +(y1-y2)**2 +(z1-z2)**2)))),\
            t7 == m.acos((x2*(x2-x3))+ (y2*(y2-y3))+(z2*(z2-z3))/ \
                  ((m.sqrt((x2**2) +(y2**2) + (z2**2)))*(m.sqrt((x2-x3)**2 \
                  +(y2-y3)**2 +(z2-z3)**2)))),\
            t8 == m.acos((x3*(x2-x3))+ (y3*(y2-y3))+(z3*(z2-z3))/ \
                  ((m.sqrt((x3**2) +(y3**2) + (z3**2)))*(m.sqrt((x2-x3)**2 \
                  +(y2-y3)**2 +(z2-z3)**2)))),\
            t9 == m.acos((x1*(x1-x3))+ (y1*(y1-y3))+(z1*(z1-z3))/ \
                  ((m.sqrt((x1**2) +(y1**2) + (z1**2)))*(m.sqrt((x1-x3)**2 \
                  +(y1-y3)**2 +(z1-z3)**2)))),\
            t10 == m.acos(((x1-x3) * (x2-x3) + (y1-y3) * (y2-y3) + (z1-z3) * (z2-z3)) / \
                  (m.sqrt((x1-x3)**2 + (y1-y3)**2 + (z1-z3)**2) * m.sqrt((x2-x3)**2 \
                  +(y2-y3)**2 + (z2-z3)**2))),\
            t11 == m.acos(((x1-x3)* (x1-x2) + (y1-y3)*(y1-y2) +(z1-z3)*(z1-z2))/ \
                  (m.sqrt((x1-x3)**2 + (y1-y3)**2 + (z1-z3)**2)* \
                  m.sqrt((x1-x2)**2 +(y1-y2)**2 +(z1-z2)**2))),\
            t12 == m.acos(((x1-x2) *(x2-x3)+(y1-y2)*(y2-y3) +(z1-z2)*(z2-z3))/ \
                  (m.sqrt((x1-x2)**2 +(y1-y2)**2 + (z1-z2)**2)* \
                  m.sqrt((x2-x3)**2 +(y2-y3)**2 +(z2-z3)**2)))])
m.open_folder()
m.options.SOLVER = 1
m.solve(disp = True)
print("vector a:",x1.value,y1.value,z1.value, \
      " vector b:",x2.value,y2.value,z2.value, \
      " vector c:",x3.value,y4.value,z3.value,)

m.open_folder()命令打开运行目录,您可以在其中找到该infeasibilities.txt文件。您也可以打开gk_model0.apm文件查看您的方程式。

于 2019-09-20T13:59:28.887 回答