Problem with the initial value in optimization problem

22 views

Kostas Mylonas

Sep 4, 2023, 10:49:37 AMSep 4
to Pyomo Forum
I have formulated the following problem that tries to minimize the cost of imported energy by controlling a battery:

```
def self_consumption(pv, demand, prices):

pv.index = np.arange(1, len(pv) + 1)
demand.index = np.arange(1, len(demand) + 1)
prices.index = np.arange(1, len(prices) + 1)

# Define the model
model = pyo.ConcreteModel()

# Define the set of timesteps
model.timesteps = pyo.Set(initialize = pyo.RangeSet(len(pv)), ordered = True)

# Define the inputs of the model
model.b_efficiency  = pyo.Param(initialize = ETTA)
model.b_cap = pyo.Param(initialize = BATTERY_CAPACITY)
model.b_min_soc = pyo.Param(initialize = BATTERY_SOC_MIN)
model.b_max_soc = pyo.Param(initialize = BATTERY_SOC_MAX)
model.b_charging_rate = pyo.Param(initialize = BATTERY_CHARGE_RATE)
model.Ppv = pyo.Param(model.timesteps, initialize = pv.to_dict()['value'], within = pyo.Any)
model.Pdemand = pyo.Param(model.timesteps, initialize = demand.to_dict()['value'], within = pyo.Any)
model.day_ahead_prices = pyo.Param(model.timesteps, initialize = prices.to_dict(), within = pyo.Any)

# Define the decision variables of the model
model.Pbat_ch = pyo.Var(model.timesteps, within = pyo.NonNegativeReals, bounds = (0, model.b_charging_rate))
model.Pbat_dis = pyo.Var(model.timesteps, within = pyo.NonNegativeReals, bounds = (0, model.b_charging_rate))
model.Ebat = pyo.Var(model.timesteps, within = pyo.NonNegativeReals, bounds = (model.b_min_soc * model.b_cap, model.b_max_soc * model.b_cap))
model.Pgrid = pyo.Var(model.timesteps)
model.is_charging = pyo.Var(model.timesteps, within = pyo.Binary)
# model.bat_ch_bin = pyo.Var(model.timesteps, within = pyo.Binary)
# model.bat_dis_bin = pyo.Var(model.timesteps, within = pyo.Binary)
model.AbsPgrid = pyo.Var(model.timesteps, within=pyo.NonNegativeReals)
# model.Y = pyo.Var(model.timesteps, within=pyo.NonNegativeReals)

# Define the constraints of the model
def BatEnergyRule(model, t):
if t == 1:
return model.Ebat[t] == model.b_cap/2 # battery initialization at half of the capacity (assumption)
else:
return model.Ebat[t] == model.Ebat[t-1] + (model.b_efficiency * model.Pbat_ch[t] - model.Pbat_dis[t]/model.b_efficiency)

model.cons1 = pyo.Constraint(model.timesteps, rule = BatEnergyRule)

def PowerBalanceRule(model, t):
return model.Pgrid[t] == model.Pdemand[t] - model.Ppv[t] + model.Pbat_ch[t] - model.Pbat_dis[t]

model.cons2 = pyo.Constraint(model.timesteps, rule = PowerBalanceRule)

# Charge/Discharge constraint
def charge_discharge_rule(model, t):
return model.Pbat_ch[t] <= model.is_charging[t] * model.b_charging_rate

model.charge_constraint = pyo.Constraint(model.timesteps, rule=charge_discharge_rule)

def discharge_charge_rule(model, t):
return model.Pbat_dis[t] <= (1 - model.is_charging[t]) * model.b_charging_rate

model.discharge_constraint = pyo.Constraint(model.timesteps, rule=discharge_charge_rule)

def AbsPgridRule1(model, t):
return model.AbsPgrid[t] >= model.Pgrid[t]

model.cons6 = pyo.Constraint(model.timesteps, rule=AbsPgridRule1)

def AbsPgridRule2(model, t):
return model.AbsPgrid[t] >= -model.Pgrid[t]

model.cons7 = pyo.Constraint(model.timesteps, rule=AbsPgridRule2)

# Define the objective function

def ObjRule(model):
return sum(model.day_ahead_prices[t] * model.AbsPgrid[t] for t in model.timesteps)

model.obj = pyo.Objective(rule=ObjRule, sense=pyo.minimize)

# Choose a solver
opt = pyo.SolverFactory('glpk')

# Solve the optimization problem
result = opt.solve(model, tee = True)

if result.solver.status != pyo.SolverStatus.ok:
logger.error(f'Solver failed to find a solution')
raise Exception('Solver failed to find a solution')

solution = {
"Pbat_ch": {t: model.Pbat_ch[t].value for t in model.timesteps},
"Pbat_dis": {t: model.Pbat_dis[t].value for t in model.timesteps},
"Ebat": {t: model.Ebat[t].value for t in model.timesteps},
"Pgrid": {t: model.Pgrid[t].value for t in model.timesteps},
}

return solution
```

However, by inspecting the solution I notice that the Ebat at the first timestep is not as I have defined it (namely b_cap/2 = 8.3/2 = 4.15, but 0.415). Do you know why is this happening?