如何解决长期的 Gekko 性能
在下面的代码中(全年 PV 斜率优化 - 每小时时间步长,CSV 数据下载 link),有什么方法可以加快优化的性能?我的建模效率低下吗?如果我在代码中将“days_to_think”变量设置为较低的天数(例如14天),则优化可以相对较快地完成,但是“days_to_think”变量增加到180+天,我的电脑就不行了t 找到解决方案。
快速获得解决方案对我来说很重要,因为我最终尝试做的是模拟微电网(光伏、建筑、电动汽车、电池等)中的最优控制。
我的代码如下所示。
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,2)
plt.step(m.time,slope_in_degree,'r')
plt.ylabel('slope [deg]')
plt.xlabel('Time [hr]')
plt.show()
解决方法
这个问题具有挑战性的部分原因是它是一个带有 max3
的混合整数非线性规划问题。 max3
函数用于在 0.45 的下限处裁剪 0.55+0.437*cos_theta+0.313*(cos_theta)**2
。
问题中的动态因素是限制角度变化速度的 DMAX
约束。您的问题中没有微分方程。如果没有 DMAX
约束,那么您可以分别解决每个时间段。对于 14 天期限,您的解决方案是:
使用 APOPT 求解器在大约 17 秒内求解。
Number of state variables: 3350
Number of total equations: - 2680
Number of slack variables: - 670
---------------------------------------
Degrees of freedom : 0
----------------------------------------------
Dynamic Control with APOPT Solver
----------------------------------------------
Iter: 1 I: 0 Tm: 16.77 NLPi: 76 Dpth: 0 Lvs: 0 Obj: -1.24E+06 Gap: 0.00E+00
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 16.814799999999998 sec
Objective : -1237075.78834978
Successful solution
---------------------------------------------------
如果您使用松弛变量和互补性约束重新表述问题,那么它的解决速度会快得多(2.9 秒)并且给出几乎相同的解决方案。
gamma = m.Var(0.5,lb=0.45)
slk = m.Var(lb=0); m.Minimize(1e-3*slk)
m.Equation(slk*(gamma-0.45)<1e-3)
m.Equation(gamma==0.55+0.437*cos_theta+0.313*(cos_theta)**2+slk)
它现在可以解决全年(365 天)的问题,并且需要 133 秒来解决具有 78,831 个变量和 61,313 个方程的问题。
Number of state variables: 78831
Number of total equations: - 61313
Number of slack variables: - 8759
---------------------------------------
Degrees of freedom : 8759
EXIT: Optimal Solution Found.
The solution was found.
The final value of the objective function is -4.464484997593126E+7
---------------------------------------------------
Solver : IPOPT (v3.12)
Solution time : 132.942 sec
Objective : -4.464484990790293E+7
Successful solution
---------------------------------------------------
这是完整的脚本。
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import csv
m = GEKKO(remote=False)
#(14 days are ok,365 days can't be solved)
days_to_consider = 365 # number of days to optimize slope
# Hourly time step
m.time = np.linspace(0,24*days_to_consider-1,\
24*days_to_consider)
# 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
# area * peak power * performance factor / 1000
misc = m.Param(value=area * P_pk * p_factor/1000)
# Initialize variables
# unit: radian
slope = m.MV(value=0.9225608396276507861,lb=0.0,ub=1.5708)
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.Var(0.5,lb=0.45)
slk = m.Var(lb=0); m.Minimize(1e-3*slk)
m.Equation(slk*(gamma-0.45)<1e-3)
m.Equation(gamma==0.55+0.437*cos_theta+0.313*(cos_theta)**2+slk)
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 = 3
#m.options.COLDSTART = 1
#m.solve(disp=True)
m.options.COLDSTART = 0
#m.options.TIME_SHIFT = 0
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(3,1,1)
plt.plot(m.time,PV_elec,'k')
plt.ylabel('PV Power [kW]')
plt.subplot(3,2)
plt.step(m.time,gamma,'b')
plt.plot([0,max(m.time)],[0.45,0.45],'k:')
plt.ylabel('gamma')
plt.subplot(3,3)
plt.step(m.time,slope_in_degree,'r')
plt.ylabel('slope [deg]')
plt.xlabel('Time [hr]')
plt.show()
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。