Question: I am working on optimizing a battery system(1 MWh) that, in conjunction with solar power(200 kW nameplate capacity), aims to reduce electricity costs for a commercial building. For those unfamiliar with commercial electricity pricing, here's a brief overview: charges typically include energy charges (based on total energy consumption over a month) and demand charges (based on the peak energy usage within a 15-minute interval, referred to as peak demand). These rates vary throughout the day based on Time of Use (TOU). The objective is to minimize monthly energy consumption to lower energy charges, while also ensuring enough energy is stored in the battery to mitigate sudden spikes in consumption to reduce demand charges. The battery should achieve this by charging when solar generation exceeds building consumption and discharging when consumption exceeds solar production. This should be straightforward with a basic load-following algorithm when sufficient solar and battery storage are available. I tested this approach with data that successfully optimized battery operations (shown in Figure 1). However, using convex optimization resulted in significantly poorer performance (Figure 2). The optimized solution from the convex solver increased energy consumption and worsened demand charges compared to not using a battery at all. Despite optimizing for TOU rates, the solver's output falls short of an ideal solution. I have thoroughly reviewed my code, objectives, and constraints, and they appear correct to me. My hypothesis is that the solver's algorithm might prioritize sending excess power to the grid (resulting in positive peaks), potentially in an an attempt negative peaks. Maybe that is why there is a random peak on the last data point. Figure 1: Near Ideal Battery Operation from the Load Following Algorithm  Figure 2: Battery Operation from the Convex Algorithm  Ideally, I aim to minimize both energy and demand charges, except when it's economical to store excess power for anticipated high-demand periods. Any insights or suggestions on refining this approach would be greatly appreciated. Thank you for your assistance. Convex Optimization CVXPy Code: # Import libraries needed import numpy as np import cvxpy as cp import matplotlib.pyplot as plt # One day of 15-minute load data load = [ 36, 42, 40, 42, 40, 44, 42, 42, 40, 32, 32, 32, 32, 30, 34, 30, 32, 30, 32, 32, 32, 32, 32, 32, 30, 32, 32, 34, 54, 62, 66, 66, 76, 76, 80, 78, 80, 80, 82, 78, 46, 104, 78, 76, 74, 78, 82, 88, 96, 84, 94, 92, 92, 92, 92, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 86, 86, 82, 66, 72, 56, 56, 54, 48, 48, 42, 50, 42, 46, 46, 46, 42, 42, 42, 44, 44, 36, 34, 32, 34, 32, 34, 32, 32] # One day of 15-minute solar data solar = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 6, 14, 26, 46, 66, 86, 104, 120, 138, 154, 168, 180, 190, 166, 152, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 190, 178, 164, 148, 132, 114, 96, 76, 58, 40, 22, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] # Define alpha matrix which are the TOU energy charges for one day lg = [31, 16, 25, 20, 4] # Length of each TOU period in 15 minute intervals pk = ['off', 'mid', 'on', 'mid', 'off'] # Classifcation of each TOU period alpha = np.array([]) for i in range(len(lg)): if pk[i] == 'on': mult = 0.1079 elif pk[i] == 'mid': mult = 0.0874 elif pk[i] == 'off': mult = 0.0755 alpha = np.append(alpha, (mult * np.ones(lg[i]))) # Define beta matricies which are the TOU demand charges for one day val = [[0.1709, 0, 0], [0, 0.0874, 0], [0, 0, 0.0755]] beta = {} for i in range(len(val)): beta_i = np.array([]) for j in range(len(lg)): if pk[j] == 'on': mult = val[0][i] elif pk[j] == 'mid': mult = val[1][i] elif pk[j] == 'off': mult = val[2][i] beta_i = np.append(beta_i, (mult * np.ones(lg[j]))) beta[i] = beta_i beta_ON = np.zeros((96, 96)) np.fill_diagonal(beta_ON, beta[0]) beta_MID = np.zeros((96, 96)) np.fill_diagonal(beta_MID, beta[1]) beta_OFF = np.zeros((96, 96)) np.fill_diagonal(beta_OFF, beta[2]) # Declare Parameters eta_plus=0.96 # charging efficiency eta_minus=0.96 # discharging efficiency Emax=900 # SOC upper limit Emin=200 # SOC lower limit E_init=500 # initial state of charge P_B_plus_max=200 # charging power limit P_B_minus_max=200 # discharging power limit opt_load=load #declaring optimal load n=96 #declaring number of timestpes for each optimization del_t=1/4 #time delta d = 1 # int(len(load) / n ) # number of days # Declare the arrays for the data outputs pg = np.array([]) psl = np.array([]) eb = np.array([]) pbp = np.array([]) pbn = np.array([]) for i in range(d): # Declare constraints List cons = [] # Pull solar and load data for nth day P_S = solar[int(n*i) : int(n*i + n)] P_L = load[int(n*i) : int(n*i + n)] # Declare variables P_G = cp.Variable(n) # Power drawn from the grid at t E_B = cp.Variable(n) # Energy in the Battery P_B_plus = cp.Variable(n) # Battery charging power at t P_B_minus = cp.Variable(n) # Battery discharging power at t P_SL = cp.Variable(n) # Solar power fed to load at t obj = cp.Minimize(cp.sum(cp.matmul(alpha, P_G) * del_t) + cp.max(cp.matmul(beta_OFF, P_G)) + cp.max(cp.matmul(beta_MID, P_G)) + cp.max(cp.matmul(beta_ON, P_G))) for t in range(n): # First iteration of constraints has an inital amount of energy for the battery. if t == 0: cons_temp = [ E_B[t] == E_init, E_B[t] >= Emin, E_B[t] <= Emax, P_B_plus[t] >= 0, P_B_plus[t] <= P_B_plus_max, P_B_minus[t] >= 0, P_B_minus[t] <= P_B_minus_max, P_SL[t] + P_B_plus[t]/eta_plus == P_S[t], P_SL[t] + P_G[t] + P_B_minus[t]*eta_minus == P_L[t], P_SL[t] >= 0 ] # Subsequent iterations use have the amount of energy from the battery calcuated from the previous constraint else: cons_temp = [ E_B[t] == E_B[t - 1] + del_t*(P_B_plus[t - 1] - P_B_minus[t - 1]), E_B[t] >= Emin, E_B[t] <= Emax, P_B_plus[t] >= 0, P_B_plus[t] <= P_B_plus_max, P_B_minus[t] >= 0, P_B_minus[t] <= P_B_minus_max, P_SL[t] + P_B_plus[t]/eta_plus == P_S[t], P_SL[t] + P_G[t] + P_B_minus[t]*eta_minus == P_L[t], P_SL[t] >= 0 ] cons += cons_temp # Solve CVX Problem prob = cp.Problem(obj, cons) prob.solve(solver=cp.CBC, verbose = True, qcp = True) # Store solution pg = np.append(pg, P_G.value) psl = np.append(psl, P_SL.value) eb = np.append(eb, E_B.value) pbp = np.append(pbp, P_B_plus.value) pbn = np.append(pbn, P_B_minus.value) # Update energy stored in battery for next iteration E_init = E_B[n - 1] # Plot Output time = np.arange(0, 24, 0.25) # 24 hours, 15-minute intervals plt.figure(figsize=(10, 6)) plt.plot(time, solar, label='Solar') plt.plot(time, [i * -1 for i in load], label='Load before Optimization') plt.plot(time, [i * -1 for i in pg], label='Load after Optimization') plt.plot(time, pbn - pbp, label='Battery Operation') # Adding labels and title plt.xlabel('Time') plt.ylabel('Demand (kW)') plt.title('Battery Optimization Output') # Adding legend plt.legend() # Display the plot plt.grid(True) plt.show()  Convex Optimization CVXPy Output: =============================================================================== CVXPY v1.3.2 =============================================================================== (CVXPY) Jun 22 03:24:36 PM: Your problem has 480 variables, 960 constraints, and 0 parameters. (CVXPY) Jun 22 03:24:36 PM: It is compliant with the following grammars: DCP, DQCP (CVXPY) Jun 22 03:24:36 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.) (CVXPY) Jun 22 03:24:36 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution. ------------------------------------------------------------------------------- Compilation ------------------------------------------------------------------------------- (CVXPY) Jun 22 03:24:36 PM: Compiling problem (target solver=CBC). (CVXPY) Jun 22 03:24:36 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> CBC (CVXPY) Jun 22 03:24:36 PM: Applying reduction Dcp2Cone (CVXPY) Jun 22 03:24:36 PM: Applying reduction CvxAttr2Constr (CVXPY) Jun 22 03:24:36 PM: Applying reduction ConeMatrixStuffing (CVXPY) Jun 22 03:24:36 PM: Applying reduction CBC (CVXPY) Jun 22 03:24:37 PM: Finished problem compilation (took 8.116e-01 seconds). ------------------------------------------------------------------------------- Numerical solver ------------------------------------------------------------------------------- (CVXPY) Jun 22 03:24:37 PM: Invoking solver CBC to obtain a solution. ------------------------------------------------------------------------------- Summary ------------------------------------------------------------------------------- (CVXPY) Jun 22 03:24:37 PM: Problem status: optimal (CVXPY) Jun 22 03:24:37 PM: Optimal value: -4.894e+01 (CVXPY) Jun 22 03:24:37 PM: Compilation took 8.116e-01 seconds (CVXPY) Jun 22 03:24:37 PM: Solver (including time spent in interface) took 5.628e-03 seconds  Load Following Algorithm Code: # Import libraries needed import matplotlib.pyplot as plt # One day of 15-minute load data load = [ 36, 42, 40, 42, 40, 44, 42, 42, 40, 32, 32, 32, 32, 30, 34, 30, 32, 30, 32, 32, 32, 32, 32, 32, 30, 32, 32, 34, 54, 62, 66, 66, 76, 76, 80, 78, 80, 80, 82, 78, 46, 104, 78, 76, 74, 78, 82, 88, 96, 84, 94, 92, 92, 92, 92, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 86, 86, 82, 66, 72, 56, 56, 54, 48, 48, 42, 50, 42, 46, 46, 46, 42, 42, 42, 44, 44, 36, 34, 32, 34, 32, 34, 32, 32] # One day of 15-minute solar data solar = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 6, 14, 26, 46, 66, 86, 104, 120, 138, 154, 168, 180, 190, 166, 152, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 190, 178, 164, 148, 132, 114, 96, 76, 58, 40, 22, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] battery = 500 output = [] # soc = [] # State of Charge of Battery net_load = [] # "Optimized Load" for i in range(96): # With non fully charged battery and excess solar: Pull power from the solar panels to charge the batteries if (battery < 900) and ((solar[i] - load[i]) >= 0): # Battery can only charge up to 100 kW if (solar[i] - load[i]) > 200: output.append(-200) else: output.append(load[i] - solar[i]) # With non depleted charged battery and excessive load: Discharge the batteries and send power to the gtid elif (battery > (200 + (load[i]/4))) and ((solar[i] - load[i]) < 0): # Battery can only discharge up to 100 kW if (solar[i] - load[i]) < -200: output.append(200) else: output.append(load[i] - solar[i]) else: output.append(0) battery += (-0.25 * output[i]) soc.append(battery / 1000) net_load.append(solar[i] - load[i] + output[i]) # Plot Output time = np.arange(0, 24, 0.25) # 24 hours, 15-minute intervals plt.figure(figsize=(10, 6)) plt.plot(time, solar, label='Solar') plt.plot(time, [i * -1 for i in load], label='Load before Optimization') plt.plot(time, net_load, label='Load after Optimization') plt.plot(time, output, label='Battery Operation') # Adding labels and title plt.xlabel('Time') plt.ylabel('Demand (kW)') plt.title('Battery Optimization Output') # Adding legend plt.legend() # Display the plot plt.grid(True) plt.show() 

Answer: A couple caveats:  I didn't look at the load following stuff, but you should compute the costs for that and compare with the optimization I think the way you are constructing the cost vector with all of the numpy/diagonals/ones/zeros/etc. is pretty confusing, but that is just an aside. I plotted the costs (alpha) and it looks fine.  First thing: You are asking CVXPY solve to use CBC solver in your solve statement. CBC is NOT a nonlinear solver and you have max() functions in your objective, which is nonlinear, unless CVXPY is doing some enormously complex linearization under the hood (unlikely) you are getting a junk result by using that solver, unless you reformulate and get rid of the max() stuff. (Another caveat: I don't know what you're doing with the max stuff in the obj, but it isn't super relevant). So try just letting CVXPY choose by omitting that flag (as I show below) After that, the results look ... well ... more credible. Over to you. I plotted the cost (x1000) so it would show, added the battery state, and flipped your battery usage line such that positive == charging (makes more sense to me.) The battery charges when rates are low and dumps what it has at peak rate, so that looks credible. Your code (w/ minor tweaks mentioned): import sys # Import libraries needed import numpy as np import cvxpy as cp import matplotlib.pyplot as plt # One day of 15-minute load data load = [36, 42, 40, 42, 40, 44, 42, 42, 40, 32, 32, 32, 32, 30, 34, 30, 32, 30, 32, 32, 32, 32, 32, 32, 30, 32, 32, 34, 54, 62, 66, 66, 76, 76, 80, 78, 80, 80, 82, 78, 46, 104, 78, 76, 74, 78, 82, 88, 96, 84, 94, 92, 92, 92, 92, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 86, 86, 82, 66, 72, 56, 56, 54, 48, 48, 42, 50, 42, 46, 46, 46, 42, 42, 42, 44, 44, 36, 34, 32, 34, 32, 34, 32, 32] # One day of 15-minute solar data solar = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 6, 14, 26, 46, 66, 86, 104, 120, 138, 154, 168, 180, 190, 166, 152, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 190, 178, 164, 148, 132, 114, 96, 76, 58, 40, 22, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] # Define alpha matrix which are the TOU energy charges for one day lg = [31, 16, 25, 20, 4] # Length of each TOU period in 15 minute intervals pk = ['off', 'mid', 'on', 'mid', 'off'] # Classifcation of each TOU period alpha = np.array([]) for i in range(len(lg)): if pk[i] == 'on': mult = 0.1079 elif pk[i] == 'mid': mult = 0.0874 elif pk[i] == 'off': mult = 0.0755 alpha = np.append(alpha, (mult * np.ones(lg[i]))) # Define beta matricies which are the TOU demand charges for one day val = [[0.1709, 0, 0], [0, 0.0874, 0], [0, 0, 0.0755]] beta = {} for i in range(len(val)): beta_i = np.array([]) for j in range(len(lg)): if pk[j] == 'on': mult = val[0][i] elif pk[j] == 'mid': mult = val[1][i] elif pk[j] == 'off': mult = val[2][i] beta_i = np.append(beta_i, (mult * np.ones(lg[j]))) beta[i] = beta_i beta_ON = np.zeros((96, 96)) np.fill_diagonal(beta_ON, beta[0]) beta_MID = np.zeros((96, 96)) np.fill_diagonal(beta_MID, beta[1]) beta_OFF = np.zeros((96, 96)) np.fill_diagonal(beta_OFF, beta[2]) # Declare Parameters eta_plus = 0.96 # charging efficiency eta_minus = 0.96 # discharging efficiency Emax = 900 # SOC upper limit Emin = 200 # SOC lower limit E_init = 500 # initial state of charge P_B_plus_max = 200 # charging power limit P_B_minus_max = 200 # discharging power limit opt_load = load # declaring optimal load n = 96 # declaring number of timestpes for each optimization del_t = 1 / 4 # time delta d = 1 # int(len(load) / n ) # number of days # Declare the arrays for the data outputs pg = np.array([]) psl = np.array([]) eb = np.array([]) pbp = np.array([]) pbn = np.array([]) print(alpha) for i in range(d): # Declare constraints List cons = [] # Pull solar and load data for nth day P_S = solar[int(n * i): int(n * i + n)] P_L = load[int(n * i): int(n * i + n)] # Declare variables P_G = cp.Variable(n) # Power drawn from the grid at t E_B = cp.Variable(n) # Energy in the Battery P_B_plus = cp.Variable(n) # Battery charging power at t P_B_minus = cp.Variable(n) # Battery discharging power at t P_SL = cp.Variable(n) # Solar power fed to load at t obj = cp.Minimize( cp.sum(cp.matmul(alpha, P_G) * del_t) + cp.max(cp.matmul(beta_OFF, P_G)) + cp.max( cp.matmul(beta_MID, P_G)) + cp.max(cp.matmul(beta_ON, P_G))) print(obj) for t in range(n): # First iteration of constraints has an inital amount of energy for the battery. if t == 0: cons_temp = [ E_B[t] == E_init, E_B[t] >= Emin, E_B[t] <= Emax, P_B_plus[t] >= 0, P_B_plus[t] <= P_B_plus_max, P_B_minus[t] >= 0, P_B_minus[t] <= P_B_minus_max, P_SL[t] + P_B_plus[t] / eta_plus == P_S[t], P_SL[t] + P_G[t] + P_B_minus[t] * eta_minus == P_L[t], P_SL[t] >= 0 ] # Subsequent iterations use have the amount of energy from the battery calcuated from the previous constraint else: cons_temp = [ E_B[t] == E_B[t - 1] + del_t * (P_B_plus[t - 1] - P_B_minus[t - 1]), E_B[t] >= Emin, E_B[t] <= Emax, P_B_plus[t] >= 0, P_B_plus[t] <= P_B_plus_max, P_B_minus[t] >= 0, P_B_minus[t] <= P_B_minus_max, P_SL[t] + P_B_plus[t] / eta_plus == P_S[t], P_SL[t] + P_G[t] + P_B_minus[t] * eta_minus == P_L[t], P_SL[t] >= 0 ] cons += cons_temp # Solve CVX Problem prob = cp.Problem(obj, cons) prob.solve(verbose=True, qcp=True) # Store solution pg = np.append(pg, P_G.value) psl = np.append(psl, P_SL.value) eb = np.append(eb, E_B.value) pbp = np.append(pbp, P_B_plus.value) pbn = np.append(pbn, P_B_minus.value) # Update energy stored in battery for next iteration E_init = E_B[n - 1] # Plot Output time = np.arange(0, 24, 0.25) # 24 hours, 15-minute intervals plt.figure(figsize=(10, 6)) plt.plot(time, solar, label='Solar') plt.plot(time, [i * -1 for i in load], label='Load before Optimization') plt.plot(time, [i * -1 for i in pg], label='Load after Optimization') plt.plot(time, pbp - pbn, label='Battery Operation') plt.plot(time, [t*10000 for t in alpha], label='cost[x10000]') plt.plot(time, eb, label='Battery state') # Adding labels and title plt.xlabel('Time') plt.ylabel('Demand (kW)') plt.title('Battery Optimization Output') # Adding legend plt.legend() # Display the plot plt.grid(True) plt.show()