TechQA.

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

1.5k views Asked by At

I'm trying to use linear optimisation to minimise the size of solar PV and battery for an off-grid property. I have solar irradiance data and household energy consumption data - I have created a year's worth (8760 data points) of the data below.

I think that this problem can solved in a linear way, however, I am seeing some strange behaviour with PulP not operating in the most optimal manner. Maybe it could be formulated better.

The amount of solar PV generated is directly proportional to the size of the PV system (PV_size) (I have assumed 20% efficiency). The solar PV output (PV_gen) and battery discharge (Pdischarge) must always meet the household demand (load). When the PV is greater than the household load, the excess PV can be used to charge the battery (Pcharge). When the excess PV is greater than the space available in the battery, we can assume that the battery gets fully charged and then the PV is curtailed. This maximum amount of charge delivered is described by Pcharge_a.

The amount discharged (Pdischarge) must be less than the available space in the battery. The battery charge state at any time is defined by Bstate[t], the maximum charge of the battery is Bmax. We can assume the battery has 100% depth of discharge and therefore it can be discharged to 0.

The objective function is minimise the cost of the system, which I have defined as the size of PV (PV_size) multiplied by cost of PV system (let's assume 500 per m2), plus the cost of the battery (let's use 1500 per kWh of battery capacity. The objective funciton is therefore to minimise:

PV_size * 500 + Bmax * 1500

I am using the PulP package with python, here is my code so far. It will return an optimal solution, however as can be seen in line 3 of the dataframe below, it creates a huge discharge and charge which is completely unecessary. I presume this is because I have not constrained how negative the discharge (Pdischarge) can be, and similarly I have not constrained the magnitude of the excess 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()
python optimization linear-programming pulp linear-optimization
Original Q&A
1

There are 1 answers

5
AirSquid On BEST ANSWER

Hello and welcome to the site. Interesting problem.

Your problem appears to be set up correctly. You are getting the "massive discharge" because you are artificially (?) constraining the battery to be at its maximum charge in the middle of the night at time=0 (assessed by looking at the periodicity of your solar energy). This would cause the battery to be artificially too large, because it doesn't need to be at max capacity at that particular time... Optimally, it should be at peak charge when the requirement becomes greater than PV generation, right?

So, it is taking a massive dump in the middle of the night to shed this charge, which (in your model) is penalty free. See the plot below from your data truncated to 5 periods:

enter image description here

So, what can be done? First, un-constrain your battery at t=0. This will allow the model to pick an optimal initial charge. You may still see anomalous behavior because the model could still keep a higher than necessary charge there and discharge as that scenario has the same objective score. So, you may choose to put a tiny penalty on cumulative discharge to influence the model. Realize this is artificially increasing the cost of battery usage very slightly, so be judicious [check this by making the discharge penalty huge in my code below and see the difference]. Or you could just ignore this, truncate the first cycle as "warm up", etc...

Here is result with no starting constraint for battery and a tiny discharge penalty...

enter image description here

Code (yours with a couple 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()

Related Questions in PYTHON

Related Questions in OPTIMIZATION

Related Questions in LINEAR-PROGRAMMING

Related Questions in PULP

Related Questions in LINEAR-OPTIMIZATION

Popular Questions

Popular Tags

javascript python java c# php android html jquery c++ css ios sql mysql r reactjs node.js arrays c asp.net json python-3.x ruby-on-rails .net sql-server swift django angular objective-c pandas excel

Trending Questions