I have a MPC problem in CVXPY for the optimal dispatch for a thermal battery and was hoping someone might have some input on how to refine my constraints or objective.

I feel like I am overcomplicated it because it seems like a simple battery dispatch problem but the primary difference is that the thermal battery can only dispatch to meet the hot water load and not the electric load. I have modeled the home electric load, the hot water load, and the PV solar generation in hourly format for all 8760 hours of the years as inputs. The home is grid-connected, so power can flow to and from the grid, but the thermal battery can only be charged from surplus or excess solar. Solar Power can flow to meet the electric load, to charge the thermal battery, or to the grid. The thermal battery discharged only to meet the hot water load.

I have tried many different objectives and penalizing and weights to try to force the thermal battery to charge but cannot figure it out. When I setup my constraints to force it, I get an infeasible solution. I think I am overthinking it or making it too complicated.

```
# Define parameters
P_PV = LA_Load['NZE_PV_Gen'].values
P_Load = LA_Load['Total Electricity (kWh)'].values
P_DHW = LA_Load['Hot Water (kWh)'].values
# Total time steps
ts = 500 # A total of 500 hours?
Horizon = 24 # Horizon for each optimization step. So hourly time steps.
E_TB_opt = np.zeros(ts + 1) # State of charge of thermal battery
P_TB_Ch_opt = np.zeros(ts) # Charging rate of thermal battery
P_TB_Dis_opt = np.zeros(ts) # Discharging rate of thermal batery
P_Grid_From_opt = np.zeros(ts) # Power from grid
P_Grid_To_opt = np.zeros(ts) # Power to grid
# Initial SOC of the thermal battery
E_TB_opt[0] = 5
# Iterate over each time step with the defined horizon
for t in range(ts - Horizon):
P_TB_Ch = cp.Variable(Horizon)
P_TB_Dis = cp.Variable(Horizon)
P_Grid_From = cp.Variable(Horizon)
P_Grid_To = cp.Variable(Horizon)
E_TB = cp.Variable(Horizon + 1)
a = cp.Variable(Horizon, boolean=True)
b = cp.Variable(Horizon, boolean=True)
constraints = [
E_TB >= 0,
E_TB <= 10.5, # Is the SOC in absolute terms instead of relative to max capacity.
P_TB_Ch >= 0,
P_TB_Ch <= 3,
P_TB_Dis >= 0,
P_TB_Dis <= 15,
P_Grid_From >= 0, # Why is both power to and from grid positive?
P_Grid_To >= 0,
a + b <= 1, # Either charging or discharging from/to the grid
P_Grid_From <= 10 * a, # P_Grid_From is non-zero only if a = 1
P_Grid_To <= 10 * b, # P_Grid_To is non-zero only if b = 1,
E_TB[0] == E_TB_opt[t], # Initial SOC for this horizon
E_TB[1:] == E_TB[:-1] + P_TB_Ch - P_TB_Dis,
P_PV[t:t + Horizon] - P_Load[t:t + Horizon] + P_Grid_From - P_Grid_To == P_TB_Ch,
P_DHW[t:t + Horizon] == P_TB_Dis
]
# Define the objective function to maximize on-site solar utilization and minimize grid power
objective = cp.Maximize(cp.sum(P_TB_Ch) - cp.sum(P_Grid_From))
# Solve the problem
prob = cp.Problem(objective, constraints)
prob.solve(solver='GLPK_MI')
if prob.status != cp.OPTIMAL:
print(f"Optimization problem at step {t} not solved to optimality: {prob.status}")
# Store optimal values for each time step within the horizon
E_TB_opt[t + 1] = E_TB.value[1] # Only update the next SOC
P_TB_Ch_opt[t] = P_TB_Ch.value[0] # Update for the first step in the horizon
P_TB_Dis_opt[t] = P_TB_Dis.value[0]
P_Grid_From_opt[t] = P_Grid_From.value[0]
P_Grid_To_opt[t] = P_Grid_To.value[0]
```