python optimization linear-programming pulp

Linear programming with PuLP - Output of variables are none

I'm trying to write a model for a nitrogen production facility to minimize the electricity costs. The facility can produce nitrogen and inject or extract it from the storage. The injection requires some additional electricity, but the extraction process does not require any energy. I have drafted the following objective function

objective function

in which the decision variable phi(i,t) denotes the flow (in m3) for the production units (phi1 - phi3) and for the injection in and extraction from the storage (phi4 & phi5). The binary variable a was put into the equation so that only one storage application (injection or extraction) is possible per stage t. electricity consumption e is a constant for each unit in kWh/m3. P(t) denotes the electricity price.

I am currently making a first version of the model with PuLP to build onto. I've tried to linearize the product of the binary variable and the continuous variables with the big M method. However, the output of the model is just 'None' for each decision variable and I can't figure out why. Looks like it cant find a solution at all. I've probably applied the big M method incorrectly. If someone could help me out that would be very nice. It's also the first piece of code I wrote so if you have any other comments please do share.

This is the program currently:

```

# Import relevant packages
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import random
import pulp as plp

# Model Creation
opt_model = plp.LpProblem(name='N2ProductionOpt')
opt_model.sense = plp.LpMinimize

# ----Input----------------------------------------------------------------------------------------------------------------------------------
# Time
set_T = list(range(0,24,1))
# Technical input
n = 3 # machine line number for machines[ASU1, ASU2, ASU3]
set_I = list(range(1, n+1))
low_cap=42000 #lower bound production capacity ASU's
max_cap = 60000 #upper bound production capacity ASU's
max_inj = max_extr = big_M = 180000 #upper bound injection/extraction capacity 
e_cons_blend = 0.25314 #electricity consumption in kWh/m3 at prod. capacity of 180.000 m3/h
e_cons_inj = 0.31045 #electricity consumption in kWh/m3 at prod. capacity of 180.000 m3/h
e_cons_extr = 0 #electricity consumption in kWh/m3
max_storage = 36.9*10**6 #max storage capacity
min_storage = 12.3*10**6 #minimal storage capacity
    
# Nitrogen demand 
n2_demand = [121548, 121453, 121537, 121715, 119228, 118547, 118675, 115909, 108003, 103060, 100284, 99211, 99915, 103157, 102453, 
             106371, 107764, 117624, 123072, 123492, 120911, 113903, 107971, 107243]
# Electricity Prices -- DA prices 
energy_prices = [107, 105, 101, 103, 109, 138, 148, 149, 144, 135, 109, 110, 111, 113, 123, 137, 147, 163, 180, 187, 148, 139, 124, 119]

#-------------------------------------------------------------------------------------------------------------------------------------------

#----Decision variables--------------------------------------------------------------------------------------------------------------------------

# production flow of each ASU
prod_flow_ASU = {(i, t): plp.LpVariable(cat='Continuous',
                                   lowBound=low_cap, upBound=max_cap, 
                                   name="x_{0}_{1}".format(i,t)) 
             for i in set_I for t in set_T}

# production flow of injection
prod_flow_inj = {t: plp.LpVariable(cat='Continuous',
                                   lowBound=0, upBound=max_inj, 
                                   name="y_{0}".format(t)) 
             for t in set_T}

# production flow of extraction
prod_flow_extr = {t: plp.LpVariable(cat='Continuous',
                                   lowBound=0, upBound=max_extr, 
                                   name="z_{0}".format(t)) 
             for t in set_T}

# amount of nitrogen available in storage
storage_level = {t: plp.LpVariable(cat='Continuous',
                                   lowBound=min_storage, upBound=max_storage, 
                                   name="s_{0}".format(t))
                 for t in set_T}

# binary value which defines the utilization, i.e. extraction or injection, of the nitrogen storage; 
storage_application = {(t): plp.LpVariable(cat='Binary',
                                   lowBound=0, upBound=1,
                                   name="l_{0}".format(t)) 
             for t in set_T} 

injection = {t: plp.LpVariable(cat='Continuous',
                                   lowBound=0, upBound=max_extr, 
                                   name="a_{0}".format(t)) 
             for t in set_T}

extraction = {t: plp.LpVariable(cat='Continuous',
                                   lowBound=0, upBound=max_extr, 
                                   name="b_{0}".format(t)) 
             for t in set_T}

# Objective function:

objective = plp.lpSum((prod_flow_ASU[i, t] * e_cons_blend + prod_flow_inj[t] * e_cons_inj + prod_flow_extr[t]*e_cons_extr) * energy_prices[t] for i in set_I for t in set_T)
opt_model.setObjective(objective)

#----Constraints-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

# Creating the binary setup of the storage utilization with the big M method
for t in set_T:
    opt_model += injection[t] <= storage_application[t] * big_M
    opt_model += injection[t] >= 0
    opt_model += injection[t] <= prod_flow_inj[t]
    opt_model += injection[t] >= prod_flow_inj[t] - (1 - storage_application[t]) * big_M

    opt_model += extraction[t] <= (1 - storage_application[t]) * big_M
    opt_model += extraction[t] >= 0
    opt_model += extraction[t] <= prod_flow_extr[t]
    opt_model += extraction[t] >= prod_flow_extr[t] - (storage_application[t]) * big_M

# Constraint to meet production demand    
    opt_model += prod_flow_ASU[1,t] + prod_flow_ASU[2,t] + prod_flow_ASU[3,t] - prod_flow_inj[t] + prod_flow_extr[t] >= n2_demand[t]
    
# Constraints for the nitrogen storage
opt_model += storage_level[0] == 36.9*10**6

for t in set_T[1:24]:
    opt_model += storage_level[t] == storage_level[t-1] + prod_flow_inj[t] - prod_flow_extr[t]
    opt_model += storage_level[t] >= 12.3*10**6
    opt_model += storage_level[t] <= 36.9*10**6
  
opt_model.solve

for t in set_T:
     print('\nFor stage {}:'.format(t))
     print('')
     for i in set_I:
             print('ASU {} flow is: {}'.format(i, prod_flow_ASU[i, t].varValue))
     print('Injection flow is: {}'.format(prod_flow_inj[t].varValue))
     print('Extraction flow is: {}'.format(prod_flow_extr[t].varValue))

```

The output is the following:

For stage 0:

ASU 1 flow is: None ASU 2 flow is: None ASU 3 flow is: None Injection flow is: None Extraction flow is: None

I have updated the objective functions and the constraints into:

# Objective function:

objective = plp.lpSum((prod_flow_ASU[i, t] * e_cons_blend + injection[t] * e_cons_inj + extraction[t]*e_cons_extr) * energy_prices[t] for i in set_I for t in set_T)
opt_model.setObjective(objective)

#----Constraints-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

# Creating the binary setup of the storage utilization with the big M method
for t in set_T:
    opt_model += injection[t] <= injection_application[t] * big_M
    opt_model += injection[t] <= prod_flow_inj[t]
    opt_model += injection[t] >= prod_flow_inj[t] - (1 - injection_application[t]) * big_M
  

    opt_model += extraction[t] <= extraction_application[t] * big_M
    opt_model += extraction[t] <= prod_flow_extr[t]
    opt_model += extraction[t] >= prod_flow_extr[t] - (1- extraction_application[t]) * big_M
   
    opt_model += extraction_application[t] + injection_application[t] <= 1
    
for i in set_I:
    for t in set_T:
        if prod_flow_ASU[i,t] != 0:
            opt_model += prod_flow_ASU[i,t] >= 42000 
        
# Constraint to meet production demand    
for t in set_T:
    opt_model += prod_flow_ASU[1,t] + prod_flow_ASU[2,t] + prod_flow_ASU[3,t] - injection[t] + extraction[t] >= n2_demand[t]
    opt_model += prod_flow_ASU[1,t] + prod_flow_ASU[2,t] + prod_flow_ASU[3,t] - n2_demand[t] == injection[t]
    
# Constraints for the nitrogen storage
opt_model += storage_level[0] == max_storage

    
for t in set_T[1:24]:

    opt_model += storage_level[t] == storage_level[t-1] + injection[t] - extraction[t]
    opt_model += storage_level[t] >= min_storage
    opt_model += storage_level[t] <= max_storage
    opt_model += storage_level[23] >= 0.98*max_storage #lower bound of 35,055 mln m3

I have the following questions:

1) I desire to constrain the variable prod_flow_ASU[i,t] to be either 0 (=off) or between 42000 and 60000 (=on). I've tried the following:

for i in set_I:
    for t in set_T:
        if prod_flow_ASU[i,t] != 0:
             opt_model += prod_flow_ASU[i,t] >= 42000

But this doesn't work. How can I model this?

2) I've tried to apply the big M method to linearize the product of the binary variable like shown in 1 & 2 , but can't seem to get it right. I've defined two binary variables, injection_application[t] and extraction_application[t] and added the constraint injection_application[t] + extraction_application[t] <= 1, so that only one operation can be applied. But I get the following output, in which the variables are not binary. When I look at my model through opt_model.solve, these variables are labeled as integers. Why is this and how can I linearize this correctly?

Big thanks for the help.

Solution

Here's a cut at this... You are close!

A couple notes:

Note the big-M constraint(s) with ASU_on. You should be using the binary variable to force the production flow to be within the min/max, or stick it to zero. Note, this also requires the ASU flow to be minimally zero, as the constraints take care of the max/min
You had some extra inject/extract variables that are not needed.... The missing piece was the ASU_on. As long as your flow-balance equation is good, the model will not arbitrarily inject/extract as that induces cost and you are minimizing cost
your storage balance equation for the first (0-th) period was incorrect. You had that the storage was full at the end, but there was "free gas" available then because there was no active flow balance, so I changed that...see the fix.
your storage capacity was so huge (I reduced it) that you could never produce and still have enough, so just think about your initial/final conditions. Perhaps you want the model to "finish where it started" (if feasible) and set the end condition to be full again (...watch the feasibility here!)... or maybe start/finish at like 80%.
As mentioned already, always check the status before you look at results
in the future, you'll likely get better help if you provide a minimally reproducible example... In this case, maybe just a couple time periods and focus on the "one thing" you are having trouble with!

Code:

# Import relevant packages
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import random
import pulp as plp

# Model Creation
opt_model = plp.LpProblem(name='N2ProductionOpt')
opt_model.sense = plp.LpMinimize

# ----Input----------------------------------------------------------------------------------------------------------------------------------
# Time
set_T = list(range(0,24,1))
# Technical input
n = 3 # machine line number for machines[ASU1, ASU2, ASU3]
set_I = list(range(1, n+1))
low_cap=42000 #lower bound production capacity ASU's
max_cap = 60000 #upper bound production capacity ASU's
max_inj = max_extr = big_M = 180000 #upper bound injection/extraction capacity 
e_cons_blend = 0.25314 #electricity consumption in kWh/m3 at prod. capacity of 180.000 m3/h
e_cons_inj = 0.31045 #electricity consumption in kWh/m3 at prod. capacity of 180.000 m3/h
e_cons_extr = 0 #electricity consumption in kWh/m3
max_storage = 125e4 #max storage capacity
min_storage = 123e3 #minimal storage capacity
    
# Nitrogen demand 
n2_demand = [121548, 121453, 121537, 0, 0, 0, 200000, 115909, 108003, 103060, 100284, 99211, 99915, 103157, 102453, 
             106371, 107764, 117624, 123072, 123492, 120911, 113903, 107971, 107243]
# Electricity Prices -- DA prices 
energy_prices = [107, 105, 101, 103, 109, 138, 148, 149, 144, 135, 109, 110, 111, 113, 123, 137, 147, 163, 180, 187, 148, 139, 124, 119]

#-------------------------------------------------------------------------------------------------------------------------------------------

#----Decision variables--------------------------------------------------------------------------------------------------------------------------

# production flow of each ASU
produce = {(i, t): plp.LpVariable(cat='Continuous',
                                   lowBound=0, upBound=max_cap,    # <-- note LB change!!
                                   name="x_{0}_{1}".format(i,t)) 
             for i in set_I for t in set_T}
# binary variable to indicate ASU is "on" and must produce minimal output
ASU_on = {(i, t): plp.LpVariable(cat='Binary', name=f'ASU_on_{i}_{t}') for i in set_I for t in set_T}


# production flow of injection
inject = {t: plp.LpVariable(cat='Continuous',
                                   lowBound=0, upBound=max_inj, 
                                   name="y_{0}".format(t)) 
             for t in set_T}

# production flow of extraction
extract = {t: plp.LpVariable(cat='Continuous',
                                   lowBound=0, upBound=max_extr, 
                                   name="z_{0}".format(t)) 
             for t in set_T}

# amount of nitrogen available in storage
storage_level = {t: plp.LpVariable(cat='Continuous',
                                   lowBound=min_storage, upBound=max_storage, 
                                   name="s_{0}".format(t))
                 for t in set_T}

# # binary value which defines the utilization, i.e. extraction or injection, of the nitrogen storage; 
# storage_application = {(t): plp.LpVariable(cat='Binary',
#                                    lowBound=0, upBound=1,
#                                    name="l_{0}".format(t)) 
#              for t in set_T} 

# injection = {t: plp.LpVariable(cat='Continuous',
#                                    lowBound=0, upBound=max_extr, 
#                                    name="a_{0}".format(t)) 
#              for t in set_T}

# extraction = {t: plp.LpVariable(cat='Continuous',
#                                    lowBound=0, upBound=max_extr, 
#                                    name="b_{0}".format(t)) 
#              for t in set_T}

# Objective function:

objective = plp.lpSum((produce[i, t] * e_cons_blend + inject[t] * e_cons_inj + extract[t]*e_cons_extr) * energy_prices[t] for i in set_I for t in set_T)
opt_model.setObjective(objective)

#----Constraints-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

# # Creating the binary setup of the storage utilization with the big M method
# for t in set_T:
#     opt_model += injection[t] <= storage_application[t] * big_M
#     opt_model += injection[t] >= 0
#     opt_model += injection[t] <= inject[t]
#     opt_model += injection[t] >= inject[t] - (1 - storage_application[t]) * big_M

#     opt_model += extraction[t] <= (1 - storage_application[t]) * big_M
#     opt_model += extraction[t] >= 0
#     opt_model += extraction[t] <= extract[t]
#     opt_model += extraction[t] >= extract[t] - (storage_application[t]) * big_M

# Constrain the min/max production by ASU and time period
for t in set_T:
    for i in set_I:
        # constrain the minimum -- force it to produce the min, if it is on
        opt_model += produce[i, t] >= ASU_on[i, t] * low_cap
        # constrain the maximum -- force it to be "on" if it produces anything, and limit it to max
        opt_model += produce[i, t] <= ASU_on[i, t] * max_cap

# Constraint to meet production demand - GOOD!
for t in set_T:
    opt_model += produce[1,t] + produce[2,t] + produce[3,t] - inject[t] + extract[t] >= n2_demand[t]
    
# Constraints for the nitrogen storage - GOOD!

for t in set_T:
    if t == 0:
        opt_model += storage_level[t] == max_storage + inject[t] - extract[t]
    else:
        opt_model += storage_level[t] == storage_level[t-1] + inject[t] - extract[t]
    # bounds are already set on storage variable, this is unneeded
    # opt_model += storage_level[t] >= 12.3*10**2
    # opt_model += storage_level[t] <= 36.9*10**6
  
opt_model.solve()

for t in set_T:
    print('\nFor stage {}:'.format(t))
    print('')
    for i in set_I:
        if ASU_on[i, t].varValue:   # will be true if binary var == 1
            print(f'ASU {i} running')
            print('   ASU {} flow is: {}'.format(i, produce[i, t].varValue))
         
    print('Injection flow is: {}'.format(inject[t].varValue))
    print('Extraction flow is: {}'.format(extract[t].varValue))
    print(f'storage at end of period is: {storage_level[t].varValue:0.2f}')

# print(opt_model)