pythonoptimizationconstraintspyomocode-documentation

Pyomo, how to make a constraint with two decision variables from two different Set()


I have an Abstract Pyomo optimization problem where I have two decision variables that are defined in two separate Set().

These are named model.BATTERY and model.HYDROGEN.

The constraint that I want to create is such that the sum of the two decision variables does not exceed a certain amount named FR. The constraint looks like this:

def LimitPshFlex_rule(model,i):
    return (sum(model.BatteryPsh[i] for i in model.BATTERY) + sum(model.HydrogenPsh[i] for i in model.HYDROGEN)) <= FR

When I call the function so that the constraint can be activated I need to input the Set() involved and the name of the rule. It looks like this:

    model.LimitPshFlex = Constraint(model.HYDROGEN, model.BATTERY, rule=LimitPshFlex_rule)

The problem with this approach is that Pyomo does not recognize it as it accepts only one Set() as input. Here is the error that I get.

TypeError: LimitPshFlex_rule() takes 2 positional arguments but 3 were given

I can also define them as two separate constraints but it saves me some coding later on if I manage to put them into one constraint so I wanted to check if anyone knows how to solve it. Here is my full code below to give you a better perspective.

 ## CREATING MODEL, SET AND PARAM
model = AbstractModel()

# Creating sets for each decision variable
model.BATTERY = Set()
model.HYDROGEN = Set()

# Creating the different variables that will be used for the constraints
model.BatteryMinPower = Param(model.BATTERY)
model.BatteryMaxPower = Param(model.BATTERY)
model.BatteryCapacity = Param(model.BATTERY)
model.BatterySoC = Param(model.BATTERY)
model.BatteryPrated = Param(model.BATTERY)
model.HydrogenMinPower = Param(model.HYDROGEN)
model.HydrogenMaxPower = Param(model.HYDROGEN)
model.HydrogenCapacity = Param(model.HYDROGEN)
model.HydrogenSoC = Param(model.HYDROGEN)
model.HydrogenPrated = Param(model.HYDROGEN)


# model.BatteryMinPower = Param(model.BATTERY, initialize=data['MinPower1'])
# model.BatteryMaxPower = Param(model.BATTERY, initialize=data['MaxPower1'])
# model.BatteryCapacity = Param(model.BATTERY, initialize=data['Capacity1'])
# model.BatterySoC = Param(model.BATTERY, initialize=data['SoC1'])
# model.BatteryPrated = Param(model.BATTERY, initialize=data['P_rated2'])
# model.HydrogenMinPower = Param(initialize=data['MinPower2'])
# model.HydrogenMaxPower = Param(initialize=data['MaxPower2'])
# model.HydrogenCapacity = Param(initialize=data['Capacity2'])
# model.HydrogenSoC = Param(initialize=data['SoC2'])
# model.HydrogenDisPrated = Param(initialize=data['P_rated2'])


# Creating dataframes that will be used for the optimization: 
PshFlex = pd.DataFrame(index=EMS.index) # Amount of power provided at each PTU 
SoC = pd.DataFrame(index=EMS.index) # The SoC of the different assets 
SoC_change = pd.DataFrame(index=EMS.index) # The change in SoC of each assets
P_rated_EMS = pd.DataFrame(index=EMS.index) # How much the asset is already in use during the request of power 
AGR_E_prog = pd.DataFrame(index=EMS.index) # The energy program proposed by the AGR. 


# Assigning the values of the EMS to the respective variables.
# Note: Soc used for defining the parameters of the optimization and SoC_change for saving the optimization results. 

SoC['battery_soc'] = EMS['battery_soc_actual']
SoC['H2_soc'] = EMS['hydrogen_soc_actual']
SoC_change['battery_soc_new'] = SoC['battery_soc'] #This colums will be the altered SoC after the bidding has been done.
SoC_change['H2_soc_new'] = SoC['H2_soc'] #This colums will be the altered SoC after the bidding has been done.
P_rated_EMS['Battery_Prated'] = -EMS['battery_actual'] 
P_rated_EMS['H2_disPrated'] = -EMS['hydrogen_actual']    

# These are just starting value to initaliaze the model but will be changed later. 
PshFlex['batt_power'] = 0 # power used for flex bid
PshFlex['H2_power'] = 0 # power used for flex bid
FR = 0 # request of power. 


# # Setting up the paramaters. A dictionary type of strucuture is used 
# # so that the values can be continuously updated:
 
data = {None: { #Somehow the Set will not properly create if I don't insert this None. I do not know the reason why :(
    'BATTERY': {None: [1]},
    'HYDROGEN': {None: [1]},
    'BatteryMaxPower': {1: 4},
    'BatteryMinPower': {1: 4},
    'BatteryCapacity': {1: 15.4}, 
    'BatterySoC': {1: SoC.iloc[0,0]}, 
    'BatteryPrated': {1: P_rated_EMS.iloc[0,0]},
    'HydrogenMinPower': {1: 0},
    'HydrogenMaxPower': {1: 2.3}, 
    'HydrogenCapacity': {1: 1998}, 
    'HydrogenSoC': {1: SoC.iloc[0,1]}, 
    'HydrogenPrated': {1: P_rated_EMS.iloc[0,1]}
}  }


# # Decision variable
model.BatteryPsh = Var(model.BATTERY, within=Reals, initialize=0)
model.HydrogenPsh = Var(model.HYDROGEN, within=Reals, initialize=0)


# Objective Function
def objective_rule(model): # For a flexibility request (FR) the power (Psh) of the assets is used to match it.  
    return FR - (sum(model.BatteryPsh[i] for i in model.BATTERY) + sum(model.HydrogenPsh[j] for j in model.HYDROGEN))



# Constraints

def MaxPowerRated_rule_BATT(model,i): # Max rated power limit
    return  model.BatteryPsh[i]  <= model.BatteryMaxPower[i]

def MinPowerRated_rule_BATT(model,i): # Min rated power limit
    return  model.BatteryMinPower[i]  <= model.BatteryPsh[i] 

def MaxCapacityLimits_rule_BATT(model,i): # Checks that the power flex is within the limits of the storage (discharge limit).
    return model.BatteryPsh[i] <= model.BatteryCapacity[i]*model.BatterySoC[i]**4   

def MaxPowerAvailable_rule_BATT(model,i): # Checks how much power the asset is already using at that moment. 
    return model.BatteryPsh[i] <= model.BatteryMaxPower[i] - model.BatteryPrated[i]

def MinPowerRated_rule_HYDRO(model,i): # Min rated power limit
    return  model.HydrogenMinPower[i] <= model.HydrogenPsh[i] 

def MaxPowerRated_rule_HYDRO(model, i): # Max rated power limit
    return  model.HydrogenPsh[i]  <= model.HydrogenMaxPower[i]

def MaxCapacityLimits_rule_HYDRO(model, i): # Checks that the power flex is within the limits of the storage (discharge limit).
    return model.HydrogenPsh[i] <= model.HydrogenCapacity[i]*model.HydrogenSoC[i]*4   

def MaxPowerAvailable_rule_HYDRO(model, i): # Checks how much power the asset is already using at that moment. 
    return model.HydrogenPsh[i] <= model.HydrogenMaxPower[i] - model.HydrogenPrated[i]

def LimitPshFlex_rule(model,i): # Makes sure that the optimization does not provide more power than what the FR needs. (minimization wil try to reduce it as much as possible)
    return (sum(model.BatteryPsh[i] for i in model.BATTERY) + sum(model.HydrogenPsh[i] for i in model.HYDROGEN)) <= FR

# def LimitPshFlex_rule_HYDRO(model,i): # Makes sure that the optimization does not provide more power than what the FR needs. (minimization wil try to reduce it as much as possible)
#     return sum(model.HydrogenPsh[i] for i in model.HYDROGEN) <= FR

# def LimitPshFlex_rule_BATT(model,i): # Makes sure that the optimization does not provide more power than what the FR needs. (minimization wil try to reduce it as much as possible)
#     return sum(model.BatteryPsh[i] for i in model.BATTERY) <= FR


# Activating the constraints
model.MaxPowerRated_BATT = Constraint(model.BATTERY, rule=MaxPowerRated_rule_BATT)
model.MinPowerRated_BATT = Constraint(model.BATTERY, rule=MinPowerRated_rule_BATT)
model.MaxCapacityLimits_BATT = Constraint(model.BATTERY, rule=MaxCapacityLimits_rule_BATT)
model.MaxPowerAvailable_BATT = Constraint(model.BATTERY, rule=MaxPowerAvailable_rule_BATT)
model.MaxPowerRated_HYDRO = Constraint(model.HYDROGEN, rule=MaxPowerRated_rule_HYDRO)
model.MinPowerRated_HYDRO = Constraint(model.HYDROGEN, rule=MinPowerRated_rule_HYDRO)
model.MaxCapacityLimits_HYDRO = Constraint(model.HYDROGEN, rule=MaxCapacityLimits_rule_HYDRO)
model.MaxPowerAvailable_HYDRO = Constraint(model.HYDROGEN, rule=MaxPowerAvailable_rule_HYDRO)
# model.LimitPshFlex_BATT = Constraint(model.BATTERY, rule=LimitPshFlex_rule_BATT)
# model.LimitPshFlex_HYDRO = Constraint(model.HYDROGEN, rule=LimitPshFlex_rule_HYDRO)
model.LimitPshFlex = Constraint(model.HYDROGEN, model.BATTERY, rule=LimitPshFlex_rule)

model.PowerProvided = Objective(rule=objective_rule, sense=minimize)
instance = model.create_instance(data)


# For loop for each PTU is created. Within the for loop the model is run through the optimization and 
# the assets change of state is updated for the following PTU. (i.e. battery discharged = SoC decreases)

for t in range(0,len(EMS.index)-1):

    # The values of the SoC and rated power are updated as the for loop iterates. 
    data = {None: { #Somehow the Set() will not create propery if I don't insert this None. I do not know the reason why :(
        'BATTERY': {None: [1]},
        'HYDROGEN': {None: [1]},
        'BatteryMaxPower': {1: 4},
        'BatteryMinPower': {1: 0},
        'BatteryCapacity': {1: 15.4}, 
        'BatterySoC': {1: SoC.iloc[t,0]}, 
        'BatteryPrated': {1: P_rated_EMS.iloc[t,0]},
        'HydrogenMinPower': {1: 0},
        'HydrogenMaxPower': {1: 2.3}, 
        'HydrogenCapacity': {1: 1998}, 
        'HydrogenSoC': {1: SoC.iloc[t,1]}, 
        'HydrogenPrated': {1: P_rated_EMS.iloc[t,1]}
    }  }

    # The flexibility requests changes at every iteration. 
    FR = FR_PTU.loc[t,'Flex_Request']
    
    # Obejective goal is set 
    model.PowerProvided = Objective(rule=objective_rule, sense=minimize)
    
    # The model is created with the updated data
    instance = model.create_instance(data) 
    
    # Checks that the battery asset is available for usage (not charging). The H2 can work also during charging because it is formed of a compressor and fuel cell.               
    if (P_rated_EMS.iloc[t,0] < 0): #or (SoC_change.iloc[t,0] < 0.2): # If assets is being charged or low in charge then disable decision variable. 
        instance.BatteryPsh[1].fix(0)
    elif (P_rated_EMS.iloc[t,0] > 0):
        instance.BatteryPsh[1].unfix()

    # if SoC_change.iloc[t,0] < 0.2: # If battery SoC is lower than 0.2 disable decision variable
    #     instance.Psh[1].fix(0)
    # elif SoC_change.iloc[t,0] > 0.2: 
    #     instance.Psh[1].unfix()
        
       
        
#     # The optimization is runned and solved.       
    opt = SolverFactory('glpk')
    opt.solve(instance)  
    print('\n \n working... \n \n ') 
    
    for i in range(0,2): 
        if i == 0:
            PshFlex.iloc[t,i] = instance.BatteryPsh[1].value 
        if i == 1:
            PshFlex.iloc[t,i] = instance.HydrogenPsh[1].value
    

Solution

  • You shouldn't be passing an index in at all in making that constraint.

    Realize when you are using the rule-function combination to make constraints, it is just a mechanism to make a similar constraint "for each" value of the index(es) that you pass in. In this case, there is only one summation constraint desired, so you don't need to pass in an index, just access the correct sets within the expression.

    Also, if it is just a single constraint desired, you can just hammer out the single expression and make the constraint from that rather than having a function do it for you. Either way. As you can see below, C1 and C2 are equivalent constraints.

    Code:

    import pyomo.environ as pyo
    
    m = pyo.ConcreteModel()
    
    m.B = pyo.Set(initialize=['b1', 'b2'], doc='set of batteries')
    m.H = pyo.Set(initialize=['h1', 'h2', 'h3'], doc='set of hydrogen cells')
    
    m.batt_use = pyo.Var(m.B)
    m.h_use    = pyo.Var(m.H)
    
    def total_use(m):
        return sum(m.batt_use[b] for b in m.B) + sum(m.h_use[h] for h in m.H) <= 50
    m.C1 = pyo.Constraint(rule=total_use)
    
    m.C2 = pyo.Constraint(expr=sum(m.batt_use[b] for b in m.B) + sum(m.h_use[h] for h in m.H) <= 50)
    
    m.pprint()
    

    Output:

    2 Set Declarations
        B : set of batteries
            Size=1, Index=None, Ordered=Insertion
            Key  : Dimen : Domain : Size : Members
            None :     1 :    Any :    2 : {'b1', 'b2'}
        H : set of hydrogen cells
            Size=1, Index=None, Ordered=Insertion
            Key  : Dimen : Domain : Size : Members
            None :     1 :    Any :    3 : {'h1', 'h2', 'h3'}
    
    2 Var Declarations
        batt_use : Size=2, Index=B
            Key : Lower : Value : Upper : Fixed : Stale : Domain
             b1 :  None :  None :  None : False :  True :  Reals
             b2 :  None :  None :  None : False :  True :  Reals
        h_use : Size=3, Index=H
            Key : Lower : Value : Upper : Fixed : Stale : Domain
             h1 :  None :  None :  None : False :  True :  Reals
             h2 :  None :  None :  None : False :  True :  Reals
             h3 :  None :  None :  None : False :  True :  Reals
    
    2 Constraint Declarations
        C1 : Size=1, Index=None, Active=True
            Key  : Lower : Body                                                            : Upper : Active
            None :  -Inf : batt_use[b1] + batt_use[b2] + h_use[h1] + h_use[h2] + h_use[h3] :  50.0 :   True
        C2 : Size=1, Index=None, Active=True
            Key  : Lower : Body                                                            : Upper : Active
            None :  -Inf : batt_use[b1] + batt_use[b2] + h_use[h1] + h_use[h2] + h_use[h3] :  50.0 :   True
    
    6 Declarations: B H batt_use h_use C1 C2