Выпуклая функция в CVXPy, похоже, не дает оптимального решения для батареи в электрической системе

Я работаю над оптимизацией аккумуляторной системы (1 МВтч), которая в сочетании с солнечной энергией (паспортная мощность 200 кВт) направлена ​​на снижение затрат на электроэнергию для коммерческого здания. Для тех, кто не знаком с коммерческими ценами на электроэнергию, вот краткий обзор: плата обычно включает в себя плату за электроэнергию (на основе общего потребления энергии за месяц) и плату за потребление (на основе пикового потребления энергии в течение 15-минутного интервала, называемого пиковым потреблением). ). Эти тарифы меняются в течение дня в зависимости от времени использования (TOU).

Цель состоит в том, чтобы свести к минимуму ежемесячное потребление энергии для снижения затрат на электроэнергию, а также обеспечить сохранение достаточного количества энергии в аккумуляторе для смягчения внезапных скачков потребления и снижения платы за потребление. Батарея должна достигать этого путем зарядки, когда выработка солнечной энергии превышает потребление здания, и разрядки, когда потребление превышает выработку солнечной энергии. Это должно быть просто с помощью базового алгоритма отслеживания нагрузки, когда имеется достаточно солнечной энергии и аккумуляторной батареи.

Я протестировал этот подход на данных, которые успешно оптимизировали работу батареи (показано на рисунке 1). Однако использование выпуклой оптимизации привело к значительному снижению производительности (рис. 2). Оптимизированное решение выпуклого решателя увеличило потребление энергии и ухудшило расходы по энергопотреблению по сравнению с полным отсутствием использования батареи. Несмотря на оптимизацию показателей TOU, результаты решателя не соответствуют идеальному решению. Я тщательно просмотрел свой код, цели и ограничения, и они кажутся мне правильными. Моя гипотеза состоит в том, что алгоритм решателя может отдавать приоритет отправке избыточной мощности в сеть (что приводит к положительным пикам), потенциально пытаясь получить отрицательные пики. Возможно, именно поэтому в последней точке данных имеется случайный пик.

Рисунок 1: Работа батареи, близкая к идеальной, по алгоритму отслеживания нагрузки Рисунок 2. Работа от батареи по выпуклому алгоритму. В идеале я стремлюсь свести к минимуму как расходы на электроэнергию, так и расходы на потребление, за исключением тех случаев, когда экономически выгодно хранить избыточную электроэнергию на ожидаемые периоды высокого спроса. Будем очень признательны за любые идеи или предложения по усовершенствованию этого подхода.

Спасибо вам за вашу помощь.

Код CVXPy выпуклой оптимизации:

# 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()

Выпуклая оптимизация CVXPy:

===============================================================================
                                     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

Загрузите следующий код алгоритма:

# 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()