Problem with the initial value in optimization problem

[Kostas Mylonas]

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?