3
votes

I know this sounds really simple, but I'm writing some code for a high school class that should work. Basically, I'm solving a pretty constrained system of equations that define a pyramid. The problem is that the solver is saying one of my equations doesn't have an equality even though it seems like they all do.

I've tried pretty much everything. I've retyped the equations and doubled checked all the syntax. I've had the same error for 2 days and I am very stuck

#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,)

The error is: Exception: @error: Equation Definition Equation without an equality (=) or inequality (>,<) false

But I don't know which equation doesn't fulfill the requirements.

1

1 Answers

0
votes

The problem is the extra brackets in your equations:

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

You can fix the problem by removing the brackets that are defining a list on the right hand side of the equation.

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

Below is a modified script that still returns an infeasible solution but with a few extra things to help you diagnose the problem.

  • Name the variables so that the infeasibilities.txt file is more readable
  • Wrap the equations with \ so they are more readable
  • Check all equations, such as the second equation, part 3. Should it be + 0.5 * m.sqrt(((x2*y3-x3*y2)**2)? It currently looks like a copy of part 2 of that equation.
  • Put realistic bounds on the variables with lb (lower bound) and ub (upper bound). For example, if it is a distance then a lower bound may be 0.01 or some other small number.
  • Include better initial guess values. Starting at zero may mean that the equations are initially evaluating at infinity or 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,)

The m.open_folder() command opens the run directory where you can find the infeasibilities.txt file. You can also open the gk_model0.apm file to see your equations.