2

I'm trying to solve non-linear equations in GEKKO using bounds and initial values. The following example illustrates two main problems I'm trying to solve for a larger system of non-linear equations.

from gekko import GEKKO
import numpy as np

m = GEKKO()

a = m.Const(value=1)
b = m.Const(value=1)
x = m.Var(value=3,lb=0)
y = m.Var(value=3,lb=0)

m.Equation(10**-x + y == a)
m.Equation(x - 10**y == b)

m.solve(disp=False)

print(x.value,y.value)

This gives a solution of:

[10.999999998] [0.99999999999]

Problem 1 - I want to solve x,y using an input range of a constants, and also not be concerned with rounding errors. For example, can I incorporate an np.linspace?

a = np.linspace(1,10,10)

Problem 2 - By the time a=4, the initial value of value=3 will not work. To get around this, I'd like to use the x,y solutions for to be the initial values for the next iteration. For example, for a=1 the solution is [10.999999998] [0.99999999999] which I'd like to be the initial values for the next iteration at a=2:

x = m.Var(value=10.999999998,lb=0)
y = m.Var(value=0.99999999999,lb=0)

and so on until it gets to a=10.

John Hedengren
  • 12,068
  • 1
  • 21
  • 25
itedin
  • 23
  • 3

1 Answers1

1

Use m.Param() for a so that a.value can be adjusted with every loop. The x.value and y.value can also be adjusted to give an initial guess. Gekko automatically stores the prior solution as the new initial guess so it is not required to reinitialize x and y.

from gekko import GEKKO
import numpy as np

m = GEKKO()

av = np.linspace(1,10,10)
xv = 3; yv = 3

a = m.Param()
b = m.Const(value=1)
x = m.Var(value=xv,lb=0)
y = m.Var(value=yv,lb=0)

m.Equation(10**-x + y == a)
m.Equation(x - 10**y == b)

for ai in av:
    a.value = ai
    x.value = xv
    y.value = yv

    m.solve(disp=False)

    # print solution
    xv = np.round(x.value[0],2)
    yv = np.round(y.value[0],2)
    print(f'a: {np.round(ai,2)} x: {xv} y: {yv}')

This gives solutions:

a: 1.0 x: 11.0 y: 1.0
a: 2.0 x: 101.0 y: 2.0
a: 3.0 x: 1001.0 y: 3.0
a: 4.0 x: 10001.0 y: 4.0
a: 5.0 x: 100001.0 y: 5.0
a: 6.0 x: 1000001.0 y: 6.0
a: 7.0 x: 10000001.0 y: 7.0
a: 8.0 x: 100000001.0 y: 8.0
a: 9.0 x: 1000000001.0 y: 9.0
a: 10.0 x: 10000000001.0 y: 10.0

Another way to solve this is with IMODE=2 to solve all of the permutations in a at once. For this problem it takes an initialization with IMODE=3. This is potentially much faster so that one larger optimization problem is solved instead of a series of optimization problems.

from gekko import GEKKO
import numpy as np

m = GEKKO()

a = m.Param(value=1)
b = m.Const(value=1)
x = m.Var(value=3,lb=0)
y = m.Var(value=3,lb=0)

m.Equation(10**-x + y == a)
m.Equation(x - 10**y == b)

m.options.IMODE=3
m.solve(disp=False)

m.options.IMODE=2
a.value=np.linspace(1,10,10)
m.solve(disp=False)

for i,ai in enumerate(a.value):
    xv = np.round(x.value[i],2)
    yv = np.round(y.value[i],2)
    print(f'a: {np.round(ai,2)} x: {xv} y: {yv}')

The result is the same:

a: 1.0 x: 11.0 y: 1.0
a: 2.0 x: 101.0 y: 2.0
a: 3.0 x: 1001.0 y: 3.0
a: 4.0 x: 10001.0 y: 4.0
a: 5.0 x: 100001.0 y: 5.0
a: 6.0 x: 1000001.0 y: 6.0
a: 7.0 x: 10000001.0 y: 7.0
a: 8.0 x: 100000001.0 y: 8.0
a: 9.0 x: 1000000001.0 y: 9.0
a: 10.0 x: 10000000001.0 y: 10.0
John Hedengren
  • 12,068
  • 1
  • 21
  • 25