In the code below (PV slope optimization over for the entire year - hourly time step, CSV data download link), is there any way to speed up the optimization's performance? Did I model inefficiently? If I set up the "days_to_consider" variable in the code as a low number of days (e.g., 14 days), the optimization can be done relatively fast, but the "days_to_consider" variable is increased to 180+ days, my computer doesn't find the solution.
It is important for me to obtain the solution fast because what I try to do eventually is to simulate the optimal control in microgrids (PV, building, EV, electric battery, etc.).
My code is shown below.
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import csv
m = GEKKO(remote=False)
days_to_consider = 14 # number of days to optimize slope (14 days are ok, 365 days can't be solved)
m.time = np.linspace(0, 24*days_to_consider-1, 24*days_to_consider) # Hourly time step
# Read the weather data from CSV
with open("PV_Input.csv", encoding='utf-8-sig') as csv_file:
csv_reader = csv.reader(csv_file, delimiter=',')
inputs = [row for row in csv.reader(csv_file)]
#initialize variables
Eb_raw = []
beta_raw = [] # unit: radians
phi_raw = []
Ed_raw = []
Et_r_raw = []
for i in range(2,24*days_to_consider+2):
Eb_raw.append(float(inputs[i][1])) # solar beam radiation
beta_raw.append(float(inputs[i][6])) # solar altitude
phi_raw.append(float(inputs[i][7])) # solar azimuth
Ed_raw.append(float(inputs[i][4])) # solar diffuse radiation
Et_r_raw.append(float(inputs[i][5])) # solar ground reflected
# Assign the time-dependent coefficients
Eb = m.Param(value=Eb_raw)
beta = m.Param(value=beta_raw)
phi = m.Param(value=phi_raw)
Ed = m.Param(value=Ed_raw)
Et_r = m.Param(value=Et_r_raw)
azimuth = 0.0 # assumed azimuth is fixed, unit: rad
rho_g = 0.14 # reflectance
area = 100 # PV area
P_pk = 250 # peak power
p_factor = 0.8 # performance factor
misc = m.Param(value=area * P_pk * p_factor/1000) # area * peak power * performance factor / 1000
# Initialize variables
slope = m.MV(value=0.9225608396276507861, lb=0.0, ub=1.5708) # unit: radian
slope.STATUS = 1
slope.DCOST = 1 # penalty for unnecessary changes
slope.DMAX = 5 # for smooth slope changes
PV_elec = m.SV()
# build PV Equation
cos_theta = m.Intermediate(m.cos(beta)*(m.cos(phi)*m.cos(azimuth)+m.sin(phi)*m.sin(azimuth))*m.sin(slope)+m.sin(beta)*m.cos(slope))
gamma = m.Intermediate(m.max3(0.45, 0.55+0.437*cos_theta+0.313*(cos_theta)**2))
m.Equation(PV_elec == misc*(Eb*(m.cos(beta)*m.cos(phi)*m.cos(azimuth)*m.sin(slope) \
+ m.cos(beta)*m.sin(phi)*m.sin(azimuth)*m.sin(slope)\
+ m.sin(beta)*m.cos(slope))\
+ Ed*(gamma*m.sin(slope) + m.cos(slope))\
+ 0.5*rho_g*(1-m.cos(slope))*(Eb*m.sin(beta)+Ed)))
m.Maximize(PV_elec)
m.options.IMODE = 6 # Optimal control
m.options.SOLVER = 1
m.solve(disp=True)
# Unit conversion to degree
conversion_rad_to_deg = 180/3.14159265359
slope_in_degree = [i*conversion_rad_to_deg for i in slope]
# Plot the results
plt.subplot(2,1,1)
plt.plot(m.time, PV_elec, 'k')
plt.ylabel('PV Power [kW]')
plt.subplot(2,1,2)
plt.step(m.time, slope_in_degree,'r')
plt.ylabel('slope [deg]')
plt.xlabel('Time [hr]')
plt.show()
