How to create a binary on-off switch for 2 variables using glpk linear programming?

Question

I have created a contrived example of my problem in an example here. It involves 2 water pipes which lead from a dam to a town. I want only 1 pipe to carry water over each timestamp, t.

I have created a linear optimisation program in Pyomo to determine how best to transfer water down the pipes, but it is failing to work based on one particular constraint - my constraint that only one pipe can transfer water at a time:

model.pipe_output1 = Var(model.T, domain = NonNegativeReals)
model.pipe_output2 = Var(model.T, domain = NonNegativeReals)


def bin_constraint1(model, t):
    return model.pipe_output2[t] == model.pipe_output2[t] - model.pipe_output1[t]
model.bin_constraint1 = Constraint(model.T, rule = bin_constraint1)

The output is 0 for both pipes in every timestamp in this model, even though the logic makes complete sense to me. What should happen is for the water to be travelling down 1 pipe only and not the other in each timestamp.

I thought my constraint makes sense because both variables are set to be non negative reals, so if the RHS is negative, then it should be floored at 0 anyway. So it forces one variable to take a value, and the other to go to 0, i.e.:

Goal result:

     pipe_output1 = [10,20,10,40,30]
     pipe_output2 = [0,0,0,0,0]

What I'm getting:

     pipe_output1 = [0,0,0,0,0]
     pipe_output2 = [0,0,0,0,0]

As a test, I used a friend's Gurobi solver license, using the constraint below, and it worked!

def bin_constraint1(model, t):
    return model.pipe_output2[t] * model.pipe_output1[t] == 0
model.bin_constraint1 = Constraint(model.T, rule = bin_constraint1)

I can't use this constraint as I don't personally have a Gurobi license. Is anyone able to help me to understand what I might be doing wrong - why does my logic fall through? Any workaround solutions?

AirSquid · Accepted Answer · 2024-07-10 15:05:40Z

The concept you are looking for here is how to enforce an either-or condition in an optimization model. So you tried 2 things...

Multiplying the variables together = 0. This works, but it is a non-linear operation, so you are giving up a lot by making your model non-linear, when this can easily be linearized
Your orig constraint. The algebra doesn't work. You essentially have:

a = a - b

which simplifies to b=0

One of the common ways to implement either-or in linear fashion is to use an indicator variable and big-M constraint. If these concepts are new, you may consider looking up some online tutorials or investing a few bucks in a used LP textbook.

This works, shows the use of big-M and switches flow based on price. Note the use of the boolean logic in C2 and C3

Code:

import pyomo.environ as pyo

### DATA

demand = {0: 15, 1: 20, 2:10}


M = max(demand.values())  # big-M value

p1_cost = {0: 2.0, 1: 1.0, 2: 5.0}
p2_cost = {0: 1.5, 1: 2.5, 2: 6.0}

### MODEL

m = pyo.ConcreteModel()

### SETS
m.T = pyo.Set(initialize=demand.keys())

### VARS
m.p1 = pyo.Var(m.T, domain=pyo.NonNegativeReals)
m.p2 = pyo.Var(m.T, domain=pyo.NonNegativeReals)

m.use_p1 = pyo.Var(m.T, domain=pyo.Binary)  # 1 if p1 in use, 0 if not

### PARAMS

m.demand = pyo.Param(m.T, initialize=demand)
m.p1_cost = pyo.Param(m.T, initialize=p1_cost)
m.p2_cost = pyo.Param(m.T, initialize=p2_cost)

### OBJ : minimize cost
m.obj = pyo.Objective(expr=sum(m.p1[t] * m.p1_cost[t] + m.p2[t] * m.p2_cost[t] for t in m.T))

### CONSTRAINTS

# meet demand
def demand(m, t):
    return m.p1[t] + m.p2[t] >= m.demand[t]
m.C1 = pyo.Constraint(m.T, rule=demand)

# only 1 pipe or other
def switch_1(m, t):
    return m.p1[t] <= m.use_p1[t] * M
m.C2 = pyo.Constraint(m.T, rule=switch_1)

def switch_2(m, t):
    return m.p2[t] <= (1 - m.use_p1[t]) * M
m.C3 = pyo.Constraint(m.T, rule=switch_2)

### CHECK IT

m.pprint()

### SOLVE IT

opt = pyo.SolverFactory('cbc')

res = opt.solve(m)

if pyo.check_optimal_termination(res):
    m.display()

Output:

1 Set Declarations
    T : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    3 : {0, 1, 2}

3 Param Declarations
    demand : Size=3, Index=T, Domain=Any, Default=None, Mutable=False
        Key : Value
          0 :    15
          1 :    20
          2 :    10
    p1_cost : Size=3, Index=T, Domain=Any, Default=None, Mutable=False
        Key : Value
          0 :   2.0
          1 :   1.0
          2 :   5.0
    p2_cost : Size=3, Index=T, Domain=Any, Default=None, Mutable=False
        Key : Value
          0 :   1.5
          1 :   2.5
          2 :   6.0

3 Var Declarations
    p1 : Size=3, Index=T
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  None :  None : False :  True : NonNegativeReals
          1 :     0 :  None :  None : False :  True : NonNegativeReals
          2 :     0 :  None :  None : False :  True : NonNegativeReals
    p2 : Size=3, Index=T
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  None :  None : False :  True : NonNegativeReals
          1 :     0 :  None :  None : False :  True : NonNegativeReals
          2 :     0 :  None :  None : False :  True : NonNegativeReals
    use_p1 : Size=3, Index=T
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  None :     1 : False :  True : Binary
          1 :     0 :  None :     1 : False :  True : Binary
          2 :     0 :  None :     1 : False :  True : Binary

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : 2.0*p1[0] + 1.5*p2[0] + p1[1] + 2.5*p2[1] + 5.0*p1[2] + 6.0*p2[2]

3 Constraint Declarations
    C1 : Size=3, Index=T, Active=True
        Key : Lower : Body          : Upper : Active
          0 :  15.0 : p1[0] + p2[0] :  +Inf :   True
          1 :  20.0 : p1[1] + p2[1] :  +Inf :   True
          2 :  10.0 : p1[2] + p2[2] :  +Inf :   True
    C2 : Size=3, Index=T, Active=True
        Key : Lower : Body                 : Upper : Active
          0 :  -Inf : p1[0] - 20*use_p1[0] :   0.0 :   True
          1 :  -Inf : p1[1] - 20*use_p1[1] :   0.0 :   True
          2 :  -Inf : p1[2] - 20*use_p1[2] :   0.0 :   True
    C3 : Size=3, Index=T, Active=True
        Key : Lower : Body                       : Upper : Active
          0 :  -Inf : p2[0] - (1 - use_p1[0])*20 :   0.0 :   True
          1 :  -Inf : p2[1] - (1 - use_p1[1])*20 :   0.0 :   True
          2 :  -Inf : p2[2] - (1 - use_p1[2])*20 :   0.0 :   True

11 Declarations: T p1 p2 use_p1 demand p1_cost p2_cost obj C1 C2 C3
Model unknown

  Variables:
    p1 : Size=3, Index=T
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :   0.0 :  None : False : False : NonNegativeReals
          1 :     0 :  20.0 :  None : False : False : NonNegativeReals
          2 :     0 :  10.0 :  None : False : False : NonNegativeReals
    p2 : Size=3, Index=T
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  15.0 :  None : False : False : NonNegativeReals
          1 :     0 :   0.0 :  None : False : False : NonNegativeReals
          2 :     0 :   0.0 :  None : False : False : NonNegativeReals
    use_p1 : Size=3, Index=T
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :   0.0 :     1 : False : False : Binary
          1 :     0 :   1.0 :     1 : False : False : Binary
          2 :     0 :   1.0 :     1 : False : False : Binary

  Objectives:
    obj : Size=1, Index=None, Active=True
        Key  : Active : Value
        None :   True :  92.5

  Constraints:
    C1 : Size=3
        Key : Lower : Body : Upper
          0 :  15.0 : 15.0 :  None
          1 :  20.0 : 20.0 :  None
          2 :  10.0 : 10.0 :  None
    C2 : Size=3
        Key : Lower : Body  : Upper
          0 :  None :   0.0 :   0.0
          1 :  None :   0.0 :   0.0
          2 :  None : -10.0 :   0.0
    C3 : Size=3
        Key : Lower : Body : Upper
          0 :  None : -5.0 :   0.0
          1 :  None :  0.0 :   0.0
          2 :  None :  0.0 :   0.0

Thanks for your awesome answer @AirSquid. I will try this now on my computer and report back. I had tried Big M with glpk solver but had no luck. Around the a = a - b logic, given that both a and b were set to non negative reals, I thought it should mean that, for example, a could equal 10 and b equal 0, or vice versa. For example, 0 - 10 = 0 works because the value of a is floored to 0 so the LHS becomes 0. Anyway, still obviously some gaps in my logic. Thanks again

Collectives™ on Stack Overflow

How to create a binary on-off switch for 2 variables using glpk linear programming?

1 Answer 1

Code:

Output:

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Code:

Output:

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related