How to use if2/3 in Gekko

Question

The problem I am optimizing is the building of power plants in a transmission network. To do this I'm placing power plants at every bus and let the optimization tell me which ones should be build to minimize running cost.

To model the placing of the plant I tried using an array of binary variables that would flag i.e. be one if the plant is used at all and 0 otherwise. Then in the Objective function to minimize I multiply this array by a constant: USEW.

I have made several attempt without any working. The one that seemed to work was using the if2 Gekko function directly in the Obj. func. However I'm getting really odd results. My code is a bit long so I'll post just the relevant lines hopefully the idea would be clear, if not please let me know and I post the whole thing.

bus=node=24
t=24
Sbase=100.
Gen = 12
VOLL = 10000.
VOLW = 50.
USEW = 100.
Pw = m.Array(m.Var,(bus,t), lb=0., ub=0., value=0.)
for b in range(bus):
    m.Minimize( np.sum(VOLL*lsh[b,:] + VOLW*Pc[b,:])*Sbase \
               + m.if2(-1.*Sbase*m.sum(Pw[b,:]),1,0)*USEW )

The problem is in the if2 part. If I remove it I get expected results but then the decision of which plant to place is lost. I tried as well with if3 but didn't work either. From what I see it seems like the optimizer is trying to minimize Pw[b,:] because the result contains only zeros. Somehow bypassing the if2 part and getting into the inner sum.

Based on the documentation, this part: m.if2(-1.*Sbase*m.sum(Pw[b,:]),1,0) should return 0 or 1 but it doesn't seem like is doing that. I'm multiplying by -1 because Pw is always positive and I want to detect when Pw>0.

I would like help on how to properly use the conditional function for this purpose. Thanks

EDIT1 Consider the following case:

from gekko import GEKKO
m = GEKKO(remote=False)
Sbase=100.
Pw = array([[[0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0],
    [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0],
    [0.0], [0.0], [0.0], [0.0], [0.0], [0.0]],
   [[10.0], [10.0], [10.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0],
    [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0],
    [0.0], [0.0], [0.0], [0.0], [0.0], [0.0]]], dtype=object)

# for np.sum(Pw[0,:])=0.0
print('sum of Pw[0,:]=', np.sum(Pw[0,:]))
print(m.if3(-1.*Sbase*np.sum(Pw[0,:]),1,0).value)
print(m.if3(-1.*Sbase*np.sum(Pw[0,:]),0,1).value)
# for np.sum(Pw[1,:])=30.0
print('sum of Pw[1,:]=', np.sum(Pw[1,:]))
print(m.if3(-1.*Sbase*np.sum(Pw[1,:]),1,0).value)
print(m.if3(-1.*Sbase*np.sum(Pw[1,:]),0,1).value)

The result is always the same: 0. It doesn't matter if I swap x1 and x2 or if the condition >=0 or <0:

0.0
sum of Pw[0,:]= 0.0
0 #result 1
0 #result 2
sum of Pw[1,:]= 30.0
0 #result 3
0 #result 4

Answer 1

One thing that you can try is to use a switch point that is 1e-3 (or a certain minimum used) instead of zero. When the switch point is at zero and the condition is 1e-10 then the output will be 1 because it is greater than the switch point. This is needed because Gekko uses gradient based optimizers that have a solution tolerance of 1e-6 (default) so a solution within that tolerance is acceptable.

There are a couple examples in the documentation that may also help. You may also want to look at the sign2/sign3 functions and the max2/max3 functions that may also give you the desired result.

if2 Documentation

IF conditional with complementarity constraint switch variable. The traditional method for IF statements is not continuously differentiable and can cause a gradient-based optimizer to fail to converge. The if2 method uses a binary switching variable to determine whether y=x1 (when condition<0) or y=x2 (when condition>=0):

if3 Documentation

IF conditional with a binary switch variable. The traditional method for IF statements is not continuously differentiable and can cause a gradient-based optimizer to fail to converge. The if3 method uses a binary switching variable to determine whether y=x1 (when condition<0) or y=x2 (when condition>=0).

Usage

y = m.if3(condition,x1,x2)

Inputs:

• condition: GEKKO variable, parameter, or expression
• x1 and x2: GEKKO variable, parameter, or expression

Output:

• y = x1 when condition<0
• y = x2 when condition>=0

from gekko import GEKKO
m = GEKKO(remote=False)
p = m.Param()
y = m.if3(p-4,p**2,p+1)

# solve with condition<0
p.value = 3
m.solve(disp=False)
print(y.value)

# solve with condition>=0
p.value = 5
m.solve(disp=False)
print(y.value)

There is additional information on logical conditions with gradient-based optimizers and the difference between 2 (MPCC) and 3 (binary) types.

Response to EDIT1

Because Gekko always uses a switch condition of 0, we can modify the switch condition with condition<swc and fit it back into the gekko form with condition-swc<0. From the example in my response, we can move the switch condition by swc=0.1.

swc = 0.1
y = m.if3(p-4-swc,p**2,p+1)

In your case, you could use swc=1e-3 or something a little higher to avoid solutions right at the switch condition. While if3 typically takes longer to solve, I generally get better results than if2, especially if there are competing objectives that interfere with the if2 MPCC.