2

I am using gekko to (try) to solve a system of equations. I run into an issue as a type of term in my equations is m.asin(q1.T@q2)/m.sqrt(1-(q1.T@q2)**2), where q1 and q2 are unit vectors in R3. When they are the same direction, the numerator and denominator become zero. When I have solved this system with fsolve I have used an if-statement to return 0 if q1.T@q2 == 1, as such that there is never an issue. As I understand here if-statements are bad, as we need the differential of the equations. I looked into m.if3 a bit, and tried that, but that doesn't help. For instance, it introduces slack variables, causing the degrees of freedom to become negative, and the system does not solve.

Any ideas how to overcome a 0 over 0 problem with gekko?

Petter
  • 81
  • 4

2 Answers2

0

You can use the m.abs3() function which provides a workaround by introducing an integer variable. This method is more suitable than the traditional method for abs as it is continuously differentiable and avoids causing the gradient-based optimizer to fail to converge.

from gekko import GEKKO
m = GEKKO()

q1 = m.Const(3)  # Unit vector in R3
q2 = m.Const(4)  # Unit vector in R3

# Using m.abs3() to handle zero over zero issue
asin_term = m.abs3(m.asin(q1.T @ q2) / m.sqrt(1 - q1.T @ q2))

# Add other equations and solve the model
m.Equations([asin_term])
m.solve()

Using m.abs3() might introduce additional integer variables, which could increase the complexity of your model. You should weigh the trade-offs between using this approach and other alternatives like if-statements, depending on your specific problem and requirements.

Kozydot
  • 515
  • 1
  • 6
0

Divide-by-zero is a common problem. Here is an example problem that illustrates the problem as x approaches 1.

from gekko import GEKKO
m = GEKKO(remote=False)
x = m.Param(0.99); y = m.Var(lb=0)
m.Equation(y==m.asin(x)/m.sqrt(1-x**2))
m.solve(disp=False)
print(y.value[0])

A first thing to try is to move the denominator to the other side of the equation so that the right-hand side doesn't go to infinity as the denominator approaches zero.

# move the denominator to the other side
m = GEKKO(remote=False)
x = m.Param(0.99); y = m.Var()
m.Equation(m.sqrt(1-x**2)*y==m.asin(x))
m.solve(disp=False)
print(y.value[0])

Gekko supports open-equations so that terms can be freely rearranged. If convergence is still a problem, add a slack variable to prevent the sqrt term from going negative.

m = GEKKO(remote=False)
x = m.Param(0.99); y = m.Var()
# minimize slack variable
s = m.Var(lb=0)
m.Equation(s>=1-x**2)
m.Minimize(s)
m.Equation(m.sqrt(s)*y==m.asin(x))
m.solve(disp=False)
print(y.value[0])

This has a lower bound of 0, but sometimes a lower bound of s>=1e-6 works better. The code segments above produce these solutions:

10.131733206
10.131733206
10.131693478
John Hedengren
  • 12,068
  • 1
  • 21
  • 25