MILP Pyomo Piecewise linear expression of battery performance

Question

I try to implement a piecewise linear expression for the battery performance within an energy optimization problem using pyomo. I am using the pyomo.Piecewise libary.

My fictive performance function of the battery charge process has following form: y = (0.6*x - 0.2*x**2 - 0.01) (y=charge-out and x=charge-in)

In the following figure the linear representation is plotted.

It has obviously a negative y intercept and therefore negative y-values at low x-values. Here comes my problem. In my mwe below the battery charge-in power (x) is forced to be about 0.0125 that charge-out power (y) is 0, in case the battery shall not be charged. Due to the negative y intercept: y(0.0125)=0.

Is there a way to implement such an equation (with negative intercept) in pyomo.Piecewise and make use of the simplified linear representation without binary variables?
Is it possible to define the equation only in the bounds where y>0? But then how to implement that the charge power can also get 0, when it is not operated?

I could not find any way to implement such an equation. In case one negelcts the y intercept the implementation is smooth.

I hope I could make my problem clear enough. Thx

import pyomo.environ as pyo
import random

##% Generate some random data for PV and Load
pv = [random.randint(0, 5) for _ in range(48)]
pv_dict = (dict(enumerate(pv,1)))

load_el = [random.randint(0, 8) for _ in range(48)]
load_el_dict = (dict(enumerate(load_el,1)))

#%% 
## Define model
model = pyo.ConcreteModel()

# Define timeperiod set
model.T = pyo.RangeSet(len(pv_dict))

# Define model parameters
model.pv = pyo.Param(model.T, initialize=pv_dict)
model.load_el= pyo.Param(model.T, initialize=load_el_dict)
model.grid_cost_buy = pyo.Param(model.T, initialize=0.4)

model.battery_eoCH = pyo.Param(initialize=1.0)
model.battery_eoDCH = pyo.Param(initialize=0.1)
model.battery_capacity = pyo.Param(initialize=5)
model.battery_power_max = pyo.Param(initialize=100)

# Define the variables             
model.battery_soc = pyo.Var(model.T, bounds=(model.battery_eoDCH, model.battery_eoCH))     # battery soct with end of ch/DCH levels 
model.grid_power_import = pyo.Var(model.T, domain=pyo.NonNegativeReals)        # grid import power 
model.grid_power_export = pyo.Var(model.T, domain=pyo.NonNegativeReals)        # grid export power 
model.battery_power_DCH = pyo.Var(model.T, domain=pyo.NonNegativeReals)        # battery discharging power
# PWA variables
model.battery_power_CH_in = pyo.Var(model.T, domain=pyo.NonNegativeReals)
model.battery_power_CH_out = pyo.Var(model.T, domain=pyo.NonNegativeReals)
model.battery_power_CH_in_norm = pyo.Var(model.T, domain=pyo.NonNegativeReals, bounds=(0,1))
model.battery_power_CH_out_norm = pyo.Var(model.T, domain=pyo.NonNegativeReals)

#%% Linearization of charge efficiency

# Define function for PWA
def f(model,t,x):
    # Normalized charge performance function y=P charge out and x=P charge in
    y = (0.6*x - 0.2*x**2 - 0.01)
    return y
# Define breakpoints
breakpoints = [0.0, 0.3, 0.6, 1.0]
# Create breakpoints dict with same index as variables
PW_PTS = {}
for idx in model.battery_power_CH_in_norm.index_set():
    PW_PTS[idx] = breakpoints
    
# Define the PWA function
model.battery_power_CH_PWA_func = pyo.Piecewise(model.T, 
                                                model.battery_power_CH_out_norm, #y
                                                model.battery_power_CH_in_norm, #x
                                                pw_pts=PW_PTS, 
                                                f_rule=f, 
                                                pw_constr_type='UB', 
                                                pw_repn='CC',
                                                force_pw=False)

# Change normalized values to absolute values
def battery_power_ch_in_rule(m,t):
    return (m.battery_power_CH_in[t] == model.battery_power_CH_in_norm[t] * m.battery_power_max)
model.battery_power_ch_in_rule_c = pyo.Constraint(model.T, rule=battery_power_ch_in_rule)

def battery_power_ch_out_rule(m,t):
    return (m.battery_power_CH_out[t] == model.battery_power_CH_out_norm[t] * m.battery_power_max)
model.battery_power_ch_out_rule_c = pyo.Constraint(model.T, rule=battery_power_ch_out_rule)

#%% Further battery constraints
# Battery SoC constraint
def battery_soc_rule(m, t):
    if t == m.T.first():
        return m.battery_soc[t] == ((m.battery_power_CH_out[t] - m.battery_power_DCH[t]) / model.battery_capacity)
    
    return m.battery_soc[t] == m.battery_soc[t-1] + ((m.battery_power_CH_out[t] - m.battery_power_DCH[t]) / model.battery_capacity)
model.battery_soc_c = pyo.Constraint(model.T, rule=battery_soc_rule)

# Define balanced electricity bus rule
def balanced_bus_rule(m, t):
    return (0 == (m.pv[t] - m.load_el[t] 
                 + m.battery_power_DCH[t] - m.battery_power_CH_in[t]
                 + m.grid_power_import[t] - m.grid_power_export[t]))
model.bus_c = pyo.Constraint(model.T, rule=balanced_bus_rule)

## Define the cost function
def obj_rule(m):
    return sum(m.grid_power_import[t]*m.grid_cost_buy[t] for t in m.T)
model.obj = pyo.Objective(rule=obj_rule, sense=1)


## Solve the problem
solver = pyo.SolverFactory('gurobi')
results = solver.solve(model)
print('Total operation costs:',pyo.value(model.obj))

MILP Pyomo Piecewise linear expression of battery performance

0 Answers0