Gekko, PWL function result is not as expected

user_4676310 0
Viewed 101

The values of z are calculated in my original code, here I've given them the values that were the result of such calculations they look like model.Equation(z[i] == someCalculation) in my original code too.

from gekko import GEKKO
model = GEKKO(remote=False)

z = model.Array(model.Var,(3))
model.Equation(z[0] == 1e-5)
model.Equation(z[1] == 4)
model.Equation(z[2] == 5.99999)

z2 = model.Array(model.Var,(3))
for i in range(3):
    d = 1e-5
    x_data = [-1e5, 0, 0+d, 6]
    y_data = [0, 0, 1, 1]
    model.pwl(z[i], z2[i], x_data, y_data)
    
model.solve()
print(z)
print(z2)

I expect z2 to be [[0] [1] [1]], but I get [[0.99035520756] [0.97749877939] [0.97771605357]], z remains the same as it should. There might be an easier way to do what I'm trying to, but I would also like to understand why this approach yields this result.

Edit

Adding an objective function- like model.Obj(model.sum(z2)) - seems to make it so the solver manipluated these variables too which is undesired to say the least. If the current setup is not capable of generating a solution I would like to see an error not the change of "fixed" variables

gekko python
1 Answers

I would expect z2 to be [[1],[1],[1]] because the PWL output steps up from 0 to 1 between 0 and +1e-5 (input). Here are some things to get closer to the solution:

  • Try the APOPT solver with m.options.SOLVER=1. It is an active set solver instead of the default IPOPT solver that uses a barrier method. The active set solver tends to return more precise solutions, especially where a barrier term is involved in defining the solution.
  • Use an integer variable for z2 (with the APOPT solver). This also will help you find an integer solution, if that is a constraint.
  • Use a tighter solution tolerance such as m.options.RTOL=1e-9. The default tolerance is 1e-6. You may also need to adjust m.options.OTOL if you need a more precise objective function.

I switched the equation to model.Equation(z[0] == -1e-5) to show a solution of [[0],[1],[1]]

from gekko import GEKKO
model = GEKKO(remote=False)

z = model.Array(model.Var,(3))
model.Equation(z[0] == -1e-5)
model.Equation(z[1] == 4)
model.Equation(z[2] == 5.99999)

z2 = model.Array(model.Var,(3)) #,integer=True)
for i in range(3):
    d = 1e-5
    x_data = [-1e5, 0, 0+d, 6]
    y_data = [0, 0, 1, 1]
    model.pwl(z[i], z2[i], x_data, y_data)

#model.options.RTOL = 1e-9
model.options.SOLVER = 1
model.solve()
print(z)
print(z2)

Edit

As you discovered, the PWL function uses slack variables. Including another objective can interfere with the correct PWL solution. The way to handle that is to may the relative objective of model.Obj(model.sum(z2)) have a smaller gradient than the gradient of the PWL function. This is an example where 0.1*sum(z2) is minimized and 0.01*sum(z) is maximized.

from gekko import GEKKO
model = GEKKO(remote=False)

z = model.Array(model.Var,(3),value=1)

z2 = model.Array(model.Var,(3)) #,integer=True)
for i in range(3):
    d = 1e-5
    x_data = [-1e5, 0, 0+d, 6]
    y_data = [0, 0, 1, 1]
    model.pwl(z[i], z2[i], x_data, y_data)

# first most important objective is PWL (gradient = 1)

# second most important objective is to minimize z2 
model.Minimize(0.1*model.sum(z2))

# third most important objective is to maximize z
model.Maximize(0.01*model.sum(z))

model.options.SOLVER = 1
model.solve()
print(z)
print(z2)

This correctly gives the solution:

[[0.0] [0.0] [0.0]]
[[0.0] [0.0] [0.0]]

Alternatively, you can maximize 0.1*sum(z2) and minimize 0.01*sum(z) with:

from gekko import GEKKO
model = GEKKO(remote=False)

z = model.Array(model.Var,(3),value=1)

z2 = model.Array(model.Var,(3)) #,integer=True)
for i in range(3):
    d = 1e-5
    x_data = [-1e5, 0, 0+d, 6]
    y_data = [0, 0, 1, 1]
    model.pwl(z[i], z2[i], x_data, y_data)

# first most important objective is PWL (gradient = 1)

# second most important objective is to minimize z2 
model.Maximize(0.1*model.sum(z2))

# third most important objective is to maximize z
model.Minimize(0.01*model.sum(z))

model.options.SOLVER = 1
model.solve()
print(z)
print(z2)

to give:

[[1e-05] [1e-05] [1e-05]]
[[1.0] [1.0] [1.0]]

I appreciate the warning message suggestion. However, it is generally too hard to find errors in models when the syntax is correct but there is a configuration that would lead to an undesirable solution. Hopefully your excellent question will help others with the PWL function.

user_23669410
Related