Replacing term in optimization leads to wrong results?

Hi there,

I am working on an "microgrid optimization for blackouts" problem. I have first built the optimization without any efficiency losses. Now, I want to introduce them. As a first step in-between, I have to model the charging and discharging of the battery (whose energy content is depicted by "e_bess") by a new variable p_bess. 

Here comes the error: I replace the charging operations (p_bess_charge[t]-p_bess_discharge[t]) with p_bess ( which is defined by: p_bess == p_bess_charge[t]-p_bess_discharge[t]).

The resulting storage size "ec_bess" before replacing it is around 1,000. After introducing p_bess, it is suddenly around 57. 

The size around 1,000 would be realistic. it should be completely the same, since I have just replaced a term.

Some help here would be incredible! :) 

I have marked the important parts in the code bold.

CODE, working:

p_bess_charge = m.addVars(range(n_timesteps), vtype = GRB.CONTINUOUS, lb=0, name = "p_bess_charge")
p_bess_discharge = m.addVars(range(n_timesteps), vtype = GRB.CONTINUOUS,lb=0, name = "p_bess_discharge")
p_diesel_discharge = m.addVars(range(n_timesteps), vtype = GRB.CONTINUOUS,lb=0, name = "p_diesel_discharge")
p_grid_final = m.addVars(range(n_timesteps), vtype = GRB.CONTINUOUS, lb=0, name = "p_grid_final")
e_bess = m.addVars(range(n_timesteps), vtype = GRB.CONTINUOUS,lb=0, name = "e_bess")
ec_bess = m.addVar(vtype=GRB.CONTINUOUS,lb=0, name="ec_bess")
p_pv = m.addVar(vtype=GRB.CONTINUOUS,lb=0, name="p_pv") # size of pv installation

m.setObjective( ((ec_bess*pr_bess*0.1)+ (pr_pv*p_pv*0.04) + quicksum(p_bess_discharge[t]*discharge_penalization+p_grid_final[t]*pr_grid*0.25+p_diesel_discharge[t]*pr_diesel*0.25 for t in range(n_timesteps))), GRB.MINIMIZE) #quicksum(... for t in range(start, end)) bildet die Summe über alle variablen von t=0 bis t=23

m.addConstr(e_bess[0] == 0)

for t in range(0,n_timesteps):
m.addConstr((p_last[t] + p_bess_charge[t] == p_bess_discharge[t]+ p_diesel_discharge[t] + pv_power[t]*p_pv + p_grid_final[t]),"a")
m.addConstr((e_bess[t] <= ec_bess),"f")
m.addConstr((p_bess_charge[t] <= ec_bess),"b")
m.addConstr((p_bess_discharge[t] <= ec_bess),"c")
if bos[t]>0:
m.addConstr((p_grid_final[t] == 0),"E")
if t > 0:
m.addConstr(e_bess[t] == (e_bess[t-1]+0.25*(p_bess_charge[t]-p_bess_discharge[t])),"x")

CODE NOT WORKING:

p_bess = m.addVars(range(n_timesteps), vtype = GRB.CONTINUOUS, name = "p_bess")
p_bess_charge = m.addVars(range(n_timesteps), vtype = GRB.CONTINUOUS, lb=0, name = "p_bess_charge")
p_bess_discharge = m.addVars(range(n_timesteps), vtype = GRB.CONTINUOUS,lb=0, name = "p_bess_discharge")
p_diesel_discharge = m.addVars(range(n_timesteps), vtype = GRB.CONTINUOUS,lb=0, name = "p_diesel_discharge")
p_grid_final = m.addVars(range(n_timesteps), vtype = GRB.CONTINUOUS, lb=0, name = "p_grid_final")
e_bess = m.addVars(range(n_timesteps), vtype = GRB.CONTINUOUS,lb=0, name = "e_bess")
ec_bess = m.addVar(vtype=GRB.CONTINUOUS,lb=0, name="ec_bess")
p_pv = m.addVar(vtype=GRB.CONTINUOUS,lb=0, name="p_pv") # size of pv installation

m.setObjective( ((ec_bess*pr_bess*0.1)+ (pr_pv*p_pv*0.04) + quicksum(p_bess_discharge[t]*discharge_penalization+p_grid_final[t]*pr_grid*0.25+p_diesel_discharge[t]*pr_diesel*0.25 for t in range(n_timesteps))), GRB.MINIMIZE) #quicksum(... for t in range(start, end)) bildet die Summe über alle variablen von t=0 bis t=23

m.addConstr(e_bess[0] == 0)

for t in range(0,n_timesteps):
m.addConstr((p_last[t] + p_bess_charge[t] == p_bess_discharge[t]+ p_diesel_discharge[t] + pv_power[t]*p_pv + p_grid_final[t]),"a")
m.addConstr((e_bess[t] <= ec_bess),"f")
m.addConstr((p_bess_charge[t] <= ec_bess),"b")
m.addConstr((p_bess_discharge[t] <= ec_bess),"c")
m.addConstr((p_bess[t] == p_bess_charge[t]-p_bess_discharge[t]),"uuuu")
if bos[t]>0:
m.addConstr((p_grid_final[t] == 0),"E")
if t > 0:
m.addConstr(e_bess[t] == (e_bess[t-1]+0.25*(p_bess[t])),"x")