Python 离网光伏和电池系统的 PulP 线性优化

Python PulP linear optimisation for off-grid PV and battery system

我正在尝试使用线性优化来最小化离网太阳能 PV 和电池的尺寸 属性。我有太阳辐照度数据和家庭能源消耗数据 - 我创建了以下一年的数据(8760 个数据点)。

我认为这个问题可以线性方式解决,但是,我看到一些奇怪的行为,PulP 没有以最佳方式运行。也许可以更好地制定。

太阳能光伏发电量与光伏系统的规模成正比 (PV_size)(我假设效率为 20%)。太阳能光伏输出 (PV_gen) 和电池放电 (Pdischarge) 必须始终满足家庭需求(负载)。当PV大于家庭负荷时,多余的PV可用于给电池充电(Pcharge)。当多余的 PV 大于电池可用的 space 时,我们可以假设电池充满电,然后 PV 被削减。此最大电荷量由 Pcharge_a 描述。

放电量 (Pdischarge) 必须小于电池中可用的 space。任意时刻的电池充电状态由Bstate[t]定义,电池的最大充电量为Bmax。我们可以假设电池有 100% 的放电深度,因此它可以放电到 0。

objective 函数是最小化系统成本,我将其定义为光伏规模 (PV_size) 乘以光伏系统成本(假设每平方米 500),加上电池的成本(让我们使用每千瓦时电池容量 1500 元。因此,objective 函数要最小化:

PV_size * 500 + Bmax * 1500

我正在使用带有 python 的 PulP 包,这是我目前的代码。它将 return 一个最佳解决方案,但是正如下面数据框的第 3 行所示,它会产生完全不必要的巨大放电和充电。我想这是因为我没有限制放电 (Pdischarge) 的负值,同样我也没有限制过剩 PV (Pcharge) 的大小。

dataframe of first few hours of operation

    load = np.array([0.580416667,0.539066667,0.390116667,0.232033333,
0.204533333,0.194716667,0.194633333,0.209233333,
0.247266668,0.407916668,0.537349998,0.576983332,
0.580216667,0.520566667,0.485200003,0.4197,
0.424300002,0.448333332,0.546983333,0.840733333,
1.320233332,0.856422014,0.921716667,0.720283335]*365)

solar_irrad = np.array([0,0,0,0,
0.846573268,6.670823882,22.34096457,48.40323145,
95.10129002,161.7686087,236.9894473,293.9150696,
305.3854497,294.6843366,251.7269744,182.2991627,
123.210826,73.11869927,33.55642336,9.910144956,
1.621109317,0.008980831,0,0]*365)

T = len(load)

# Decision variables
Bmax = LpVariable('Bmax', 0, None) # battery max energy (kWh)
PV_size = LpVariable('PV_size', 0, None) # PV size

# Optimisation problem
prb = LpProblem('Battery_Operation', LpMinimize)

# Objective function
prb += (PV_size*500) + (Bmax*1500)  # cost of battery

# Auxilliary variables
PV_gen = [LpVariable('PV_gen_{}'.format(i), 0, None) for i in range(T)]

# Load difference
Pflow = [LpVariable('Pflow_{}'.format(i), None, None) for i in range(T)]
# Excess PV
Pcharge = [LpVariable('Pcharge_{}'.format(i), lowBound=0, upBound=None) for i in range(T)]
# Discharge required
Pdischarge = [LpVariable('Pdischarge_{}'.format(i), lowBound=None, upBound=0) for i in range(T)]
# Charge delivered
Pcharge_a = [LpVariable('Pcharge_a{}'.format(i), 0, None) for i in range(T)]

# Battery
Bstate = [LpVariable('E_{}'.format(i), 0, None) for i in range(T)]

# Battery Constraints
prb += Bstate[0] == Bmax + Pdischarge[0] + Pcharge_a[0]
for t in range(1, T):
    prb += Bstate[t] == Bstate[t-1] + Pdischarge[t] + Pcharge_a[t] 

# Power flow Constraints
for t in range(0, T):
    
    # PV generation
    prb += PV_gen[t] == PV_size*0.2*solar_rad[t]/1000
    
    # Pflow is the energy flow reuired to meet the load
    # Negative if load greater than PV, positive if PV greater than load
    prb += Pflow[t] == PV_gen[t] - load[t]
    
    # Given the below, it will push Pflow available for charge to zero or to to or greater than excess PV
    prb += Pcharge[t] >= 0
    prb += Pcharge[t] >= Pflow[t]

    # If Pflow is negative (discharge), then it will at least ePflowual discharge rePflowuired
    # If Pflow is positive (charge), then Pdischarge (discharge rePflowuired will ePflowual 0)
    prb += Pdischarge[t] <= 0
    prb += Pdischarge[t] <= Pflow[t]
    # Discharge cannot exceed available charge in battery
    # Discharge is negative
    prb += Pdischarge[t] >= (-1)*Bstate[t-1]
    
    # Ensures that energy flow rePflowuired is satisifed by charge and discharge flows
    prb += Pflow[t] == Pcharge[t] + Pdischarge[t] 
    
    # Limit amount charge delivered by the available space in the battery
    prb += Pcharge_a[t] >= 0
    prb += Pcharge_a[t] <= Pcharge[t]
    prb += Pcharge_a[t] <= Bmax - Bstate[t-1]
    
    prb += Bstate[t] >= 0
    prb += Bstate[t] <= Bmax
    
    # Solve problem
    prb.solve()

您好,欢迎来到本网站。有趣的问题。

您的问题似乎设置正确。您正在获得“大量放电”,因为您人为地(?)将电池限制在半夜 time=0 时处于最大电量(通过查看太阳能的周期性进行评估)。这会导致电池人为地过大,因为它不需要在那个特定时间处于最大容量......最佳情况下,当需求变得大于 PV 发电时,它应该处于峰值充电,对吧?

所以,它在半夜进行大量倾倒以摆脱这种负担,这(在您的模型中)是免费的。从截断为 5 个周期的数据中查看下图:

那么,可以做什么呢?首先,un-constrain 你的电池在 t=0。这将允许模型选择最佳的初始电荷。您可能仍然会看到异常行为,因为模型 可能 仍然保持高于必要的充电和放电,因为该场景具有相同的 objective 分数。因此,您可以选择对累积排放施加 tiny 惩罚以影响模型。意识到这是人为地稍微增加了电池使用成本,所以要谨慎[通过在我下面的代码中使放电惩罚变得巨大来检查这一点并查看差异]。或者您可以忽略它,将第一个周期截断为“热身”等...

这是没有电池启动限制和微小放电惩罚的结果...

代码(你和一对 tweaks/comments):

from pulp import *
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

load = np.array([0.580416667,0.539066667,0.390116667,0.232033333,
0.204533333,0.194716667,0.194633333,0.209233333,
0.247266668,0.407916668,0.537349998,0.576983332,
0.580216667,0.520566667,0.485200003,0.4197,
0.424300002,0.448333332,0.546983333,0.840733333,
1.320233332,0.856422014,0.921716667,0.720283335]*5)

solar_rad = np.array([0,0,0,0,
0.846573268,6.670823882,22.34096457,48.40323145,
95.10129002,161.7686087,236.9894473,293.9150696,
305.3854497,294.6843366,251.7269744,182.2991627,
123.210826,73.11869927,33.55642336,9.910144956,
1.621109317,0.008980831,0,0]*5)

T = len(load)

# Decision variables
Bmax = LpVariable('Bmax', 0, None) # battery max energy (kWh)
PV_size = LpVariable('PV_size', 0, None) # PV size

# Optimisation problem
prb = LpProblem('Battery_Operation', LpMinimize)

# Auxilliary variables
PV_gen = [LpVariable('PV_gen_{}'.format(i), 0, None) for i in range(T)]

# Load difference
Pflow = [LpVariable('Pflow_{}'.format(i), None, None) for i in range(T)]
# Excess PV
Pcharge = [LpVariable('Pcharge_{}'.format(i), lowBound=0, upBound=None) for i in range(T)]
# Discharge required
Pdischarge = [LpVariable('Pdischarge_{}'.format(i), lowBound=None, upBound=0) for i in range(T)]
# Charge delivered
Pcharge_a = [LpVariable('Pcharge_a{}'.format(i), 0, None) for i in range(T)]

###  Moved this down as it needs to include Pdischarge
# Objective function
# cost + some small penalty for cumulative discharge, just to shape behavior
prb += (PV_size*500) + (Bmax*1500)  - 0.01 * lpSum(Pdischarge[t] for t in range(T)) 

# Battery
Bstate = [LpVariable('E_{}'.format(i), 0, None) for i in range(T)]

# Battery Constraints
### NOTE this is killed to allow initial state to "float"
#prb += Bstate[0] == Bmax + Pdischarge[0] + Pcharge_a[0]
for t in range(1, T):
    prb += Bstate[t] == Bstate[t-1] + Pdischarge[t] + Pcharge_a[t] 

# Power flow Constraints
for t in range(0, T):
    
    # PV generation
    prb += PV_gen[t] == PV_size*0.2*solar_rad[t]/1000
    
    # Pflow is the energy flow reuired *from the battery* to meet the load
    # Negative if load greater than PV, positive if PV greater than load
    prb += Pflow[t] == PV_gen[t] - load[t]
    
    # Given the below, it will push Pflow available for charge to zero or to to or greater than excess PV
    prb += Pcharge[t] >= 0
    prb += Pcharge[t] >= Pflow[t]

    # If Pflow is negative (discharge), then it will at least ePflowual discharge rePflowuired
    # If Pflow is positive (charge), then Pdischarge (discharge rePflowuired will ePflowual 0)
    prb += Pdischarge[t] <= 0
    prb += Pdischarge[t] <= Pflow[t]
    # Discharge cannot exceed available charge in battery
    # Discharge is negative
    prb += Pdischarge[t] >= (-1)*Bstate[t-1]
    
    # Ensures that energy flow rePflowuired is satisifed by charge and discharge flows
    prb += Pflow[t] == Pcharge[t] + Pdischarge[t] 
    
    # Limit amount charge delivered by the available space in the battery
    prb += Pcharge_a[t] >= 0
    prb += Pcharge_a[t] <= Pcharge[t]
    prb += Pcharge_a[t] <= Bmax - Bstate[t-1]
    
    prb += Bstate[t] >= 0
    prb += Bstate[t] <= Bmax
    
# Solve problem
prb.solve()

# make some records to prep for dataframe (what a pain in pulp!!)
res = []
for t in range(T):
    record = {  'period': t,
                'Load': load[t],
                'PV_gen': PV_gen[t].varValue,
                'Pflow' : Pflow[t].varValue,
                'Pcharge': Pcharge[t].varValue,
                'Pcharge_a': Pcharge_a[t].varValue,
                'Pdischarge': Pdischarge[t].varValue,
                'Bstate': Bstate[t].varValue}
    res.append(record)

df = pd.DataFrame.from_records(res)
df.set_index('period', inplace=True)
df = df.round(2)
print(df.to_string())

print(f'PV size: {PV_size.varValue : 0.1f}, Batt size: {Bmax.varValue : 0.1f}')

df.plot()
plt.show()