在下面的代码中(全年的 PV 斜率优化 - 每小时时间步长，CSV 数据下载 link )，有什么方法可以加快优化的性能吗？我建模效率低下吗？如果我在代码中将“days_to_consider”变量设置为较低的天数(例如14天)，则优化可以相对较快地完成，但是“days_to_consider”变量增加到180+天，我的计算机没有找不到解决方案。
快速获得解决方案对我来说很重要，因为我最终要做的是模拟微电网(光伏、建筑、电动汽车、电池等)中的最优控制。
我的代码如下所示。
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()
最佳答案
该问题之所以具有挑战性，部分原因在于它是一个具有 max3 的混合整数非线性规划问题。 max3 函数用于将 0.55+0.437*cos_theta+0.313*(cos_theta)**2 裁剪在下限 0.45 处。
问题中的动态是 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,1,2)
plt.step(m.time, gamma,'b')
plt.plot([0,max(m.time)],[0.45,0.45],'k:')
plt.ylabel('gamma')
plt.subplot(3,1,3)
plt.step(m.time, slope_in_degree,'r')
plt.ylabel('slope [deg]')
plt.xlabel('Time [hr]')
plt.show()
关于Gekko 长期性能，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/65903805/