[FIXED] How to use binary variable in another constraint in PYOMO?

Issue

I have created an RTC energy dispatch model using PYOMO. There is a demand profile and to serve this demand profile we have three main components: Solar, Wind, and Battery just so that the model won't give infeasible errors in case there is not enough energy supply from the three main sources I have added one backup generator called Lost_Load I have added one binary variable called insufficient_dispatch which will be 1 if lost_load >= 10% of demand profile else 0. Here is the inputs that I am using:

Here is the code:

import datetime
import pandas as pd
import numpy as np
from pyomo.environ import *




profiles = pd.read_excel('profiles.xlsx', sheet_name='15min', usecols='B:D')
solar_profile = profiles.iloc[:, 0].values
wind_profile = profiles.iloc[:, 1].values
load_profile = profiles.iloc[:, 2].values


interval_freq = 60
solar_capacity = 120
wind_capacity = 170
power = 90
energy = 360

model = ConcreteModel()

model.m_index = Set(initialize=list(range(len(load_profile))))
# model.total_instance = Param(initialize=35010000)



model.chargeEff = Param(initialize=0.92)
model.dchargeEff = Param(initialize=0.92)
model.storage_power = Param(initialize=power)
model.storage_energy = Param(initialize=energy)

#variable
# model.grid = Var(model.m_index, domain=Reals)
model.grid = Var(model.m_index, domain=NonNegativeReals)
model.p_max = Var(domain=NonNegativeReals)

# storage
#variable
model.e_in = Var(model.m_index, domain=NonNegativeReals)
model.e_out = Var(model.m_index, domain=NonNegativeReals)
model.soc = Var(model.m_index, domain=NonNegativeReals)
model.ein_eout_lmt = Var(model.m_index, domain=Binary)

#  Solar variable

solar_ac = np.minimum(solar_profile * solar_capacity * (interval_freq / 60), solar_capacity * (interval_freq / 60))
model.solar_cap = Param(initialize=solar_capacity)
model.solar_use = Var(model.m_index, domain=NonNegativeReals)

#  Wind variables

wind_ac = np.minimum(wind_profile * wind_capacity * (interval_freq / 60), wind_capacity * (interval_freq / 60))
model.wind_cap = Param(initialize=wind_capacity)
model.wind_use = Var(model.m_index, domain=NonNegativeReals)

# Solar profile
model.solar_avl = Param(model.m_index, initialize=dict(zip(model.m_index, solar_ac)))

# Wind profile
model.wind_avl = Param(model.m_index, initialize=dict(zip(model.m_index, wind_ac)))

# Load profile
model.load_profile = Param(model.m_index, initialize=dict(zip(model.m_index, load_profile * interval_freq / 60.0)))

model.lost_load = Var(model.m_index, domain=NonNegativeReals)

# Auxiliary binary variable for the condition
model.insufficient_dispatch = Var(model.m_index, domain=Binary)


# Objective function
def revenue(model):
    total_revenue = sum(
        model.grid[m] * 100 +
        model.lost_load[m] * -1000000000 #* model.insufficient_dispatch[m]
        for m in model.m_index)

    return total_revenue


model.obj = Objective(rule=revenue, sense=maximize)


eps = 1e-3  # to create a gap for gen3_status constraints
bigm = 1e3  # choose this value high but not so much to avoid numeric instability

# Create 2 BigM constraints
def gen3_on_off1(model, m):
    return model.lost_load[m] >= 0.1 * model.load_profile[m] + eps - bigm * (1 - model.insufficient_dispatch[m])

def gen3_on_off2(model, m):
    return model.lost_load[m] <= 0.1 * model.load_profile[m] + bigm * model.insufficient_dispatch[m]

model.gen3_on_off1 = Constraint(model.m_index, rule=gen3_on_off1)
model.gen3_on_off2 = Constraint(model.m_index, rule=gen3_on_off2)


def energy_balance(model, m):
    return model.grid[m] <= model.solar_use[m] + model.wind_use[m] + model.e_out[m] - model.e_in[m] + model.lost_load[m]


model.energy_balance = Constraint(model.m_index, rule=energy_balance)


def grid_limit(model, m):
    return model.grid[m] >= model.load_profile[m]

model.grid_limit = Constraint(model.m_index, rule=grid_limit)


def max_solar_gen(model, m):
    eq = model.solar_use[m] <= model.solar_avl[m] * 1
    return eq


model.max_solar_gen = Constraint(model.m_index, rule=max_solar_gen)


def min_solar_gen(model, m):
    eq = model.solar_use[m] >= model.solar_avl[m] * 0.6
    return eq


# model.min_solar_gen = Constraint(model.m_index, rule=min_solar_gen)


def max_wind_gen(model, m):
    eq = model.wind_use[m] <= model.wind_avl[m] * 1
    return eq


model.max_wind_gen = Constraint(model.m_index, rule=max_wind_gen)


def min_wind_gen(model, m):
    eq = model.wind_use[m] >= model.wind_avl[m] * 0.3
    return eq


# model.min_wind_gen = Constraint(model.m_index, rule=min_wind_gen)


# Charging and discharging controller
def ein_limit1(model, m):
    eq = model.e_in[m] + (1 - model.ein_eout_lmt[m]) * 1000000 >= 0

    return eq

model.ein_limit1 = Constraint(model.m_index, rule=ein_limit1)


def eout_limit1(model, m):
    eq = model.e_out[m] + (1 - model.ein_eout_lmt[m]) * -1000000 <= 0
    return eq

model.eout_limit1 = Constraint(model.m_index, rule=eout_limit1)


def ein_limit2(model, m):
    eq = model.e_out[m] + (model.ein_eout_lmt[m]) * 1000000 >= 0
    return eq

model.ein_limit2 = Constraint(model.m_index, rule=ein_limit2)


def eout_limit2(model, m):
    eq = model.e_in[m] + (model.ein_eout_lmt[m]) * -1000000 <= 0
    return eq

model.eout_limit2 = Constraint(model.m_index, rule=eout_limit2)


#max charging
def max_charge(model, m):
    return model.e_in[m] <= model.storage_power*interval_freq/60

model.max_charge = Constraint(model.m_index, rule=max_charge)


# max discharging
def max_discharge(model, m):
    return model.e_out[m] <= model.storage_power*interval_freq/60

model.max_max_discharge = Constraint(model.m_index, rule=max_discharge)


def soc_update(model, m):
    if m == 0:
        eq = model.soc[m] == \
             (
                     (model.storage_energy * 1) +
                     (model.e_in[m] * model.chargeEff) -
                     (model.e_out[m] / model.dchargeEff) #* model.insufficient_dispatch[m]
             )
    else:
        eq = model.soc[m] == \
             (
                     model.soc[m - 1] +
                     (model.e_in[m] * model.chargeEff) -
                     (model.e_out[m] / model.dchargeEff) #* model.insufficient_dispatch[m]
             )

    return eq


model.soc_update = Constraint(model.m_index, rule=soc_update)


def soc_min(model, m):
    eq = model.soc[m] >= model.storage_energy * 0.2

    return eq


model.soc_min = Constraint(model.m_index, rule=soc_min)


def soc_max(model, m):
    eq = model.soc[m] <= model.storage_energy * 1

    return eq


model.soc_max = Constraint(model.m_index, rule=soc_max)

def throughput_limit(model):
    energy_sum = sum(model.e_out[idx] for idx in model.m_index)
    return energy_sum <= model.storage_energy*365


# model.throughput_limit = Constraint(rule=throughput_limit)


Solver = SolverFactory('gurobi')
Solver.options['LogFile'] = "gurobiLog"
Solver.options['MIPGap'] = 0.50
print('\nConnecting to Gurobi Server...')
results = Solver.solve(model)

if (results.solver.status == SolverStatus.ok):
    if (results.solver.termination_condition == TerminationCondition.optimal):
        print("\n\n***Optimal solution found***")
        print('obj returned:', round(value(model.obj), 2))
    else:
        print("\n\n***No optimal solution found***")
        if (results.solver.termination_condition == TerminationCondition.infeasible):
            print("Infeasible solution")
            exit()
else:
    print("\n\n***Solver terminated abnormally***")
    exit()

grid_use = []
solar = []
wind = []
e_in = []
e_out=[]
soc = []
lost_load = []
insufficient_dispatch = []
load = []

for i in range(len(load_profile)):
    grid_use.append(value(model.grid[i]))
    solar.append(value(model.solar_use[i]))
    wind.append(value(model.wind_use[i]))
    lost_load.append(value(model.lost_load[i]))
    e_in.append(value(model.e_in[i]))
    e_out.append(value(model.e_out[i]))
    soc.append(value(model.soc[i]))
    load.append(value(model.load_profile[i]))
    insufficient_dispatch.append(value(model.insufficient_dispatch[i]))

df_out = pd.DataFrame()

df_out["Grid"] = grid_use
df_out["Potential Solar"] = solar_ac
df_out["Potential Wind"] = wind_ac
df_out["Actual Solar"] = solar
df_out["Actual Wind"] = wind
# df_out["Solar Curtailment"] = solar_ac - np.array(solar)
# df_out["Wind Curtailment"] = wind_ac - np.array(wind)
df_out["Charge"] = e_in
df_out["Discharge"] = e_out
df_out["soc"] = soc
df_out["Lost Load"] = lost_load
df_out["Load profile"] = load
df_out["Dispatch"] = df_out["Grid"] - df_out["Lost Load"]
df_out["Insufficient Supply"] = insufficient_dispatch

summary = {
    'Grid': [df_out['Grid'].sum()/1000],
    'Wind': [df_out['Actual Wind'].sum()/1000],
    'Solar': [df_out['Actual Solar'].sum()/1000],
    'Storage Discharge': [df_out['Discharge'].sum()/1000],
    'Storage Charge': [df_out['Charge'].sum()/1000],
    'Shortfall': [df_out['Lost Load'].sum()/1000],
    # 'Solar Curtailment': [df_out['Solar Curtailment'].sum()/1000],
    # 'Wind Curtailment': [df_out['Wind Curtailment'].sum()/1000],
}


summaryDF = pd.DataFrame(summary)

timestamp = datetime.datetime.now().strftime('%y-%m-%d_%H-%M')

with pd.ExcelWriter('output solar '+str(solar_capacity)+'MW wind '+str(wind_capacity)+'MW ESS '+str(power)+'MW and '+str(energy)+'MWh '+str(timestamp)+'.xlsx') as writer:
    df_out.to_excel(writer, sheet_name='dispatch')
    summaryDF.to_excel(writer, sheet_name='Summary', index=False)

I want to use this binary variable insufficient_dispatch with another constraint energy_balance ( which is already there in the code ) :

def energy_balance(model, m):
    return model.grid[m] <= model.solar_use[m] + model.wind_use[m] + model.e_out[m] - model.e_in[m] + model.lost_load[m]

it should be like:

if model.insufficient_dispatch[m] == 0:
    return model.grid[m] <= model.solar_use[m] + model.wind_use[m] + model.e_out[m] - model.e_in[m] + model.lost_load[m]
else:
    return model.grid[m] * 0.9 <= model.solar_use[m] + model.wind_use[m] + model.e_out[m] - model.e_in[m] + model.lost_load[m]

But this won't work, can someone please help? This is the graphical representation of one of the months:

how i want is something like this:

Solution

If we are only interested in the left hand side of your inequality, you want 1 * grid when insufficient_dispatch == 0 or 0.9 * grid otherwise.

if i == 0:
  1.0 * x <= F
else:
  0.9 * x <= F

is equivalent to:

(1 - w) + w * 0.9

so:

w=0, (1 - 0) + 0 * 0.9 gives 1.0
w=1, (1 - 1) + 1 * 0.9 gives 0.9

You can rewrite your constraint like:

def energy_balance(model, m):
    lhs = ((1 - model.insufficient_dispatch[m]) + model.insufficient_dispatch[m] * 0.9) * model.grid[m]
    return lhs <= model.solar_use[m] + model.wind_use[m] + model.e_out[m] - model.e_in[m] + model.lost_load[m]

Edit

To keep your model linear, create 2 constraints (same explanation than my answer of your previous question):

bigm = 1e3

# When i=0
def energy_balance1(model, m):
   return model.grid[m] <= model.solar_use[m] + model.wind_use[m] + model.e_out[m] - model.e_in[m] + model.lost_load[m] + model.insufficient_dispatch[m] * bigm

# When i=1
def energy_balance2(model, m):
    return model.grid[m] * 0.9 <= model.solar_use[m] + model.wind_use[m] + model.e_out[m] - model.e_in[m] + model.lost_load[m] + (1 - model.insufficient_dispatch[m]) * bigm

Answered By - Corralien