Gekko 长期性能

在下面的代码中（全年 PV 斜率优化 - 每小时时间步长，CSV 数据下载link https://drive.google.com/file/d/172HaD87t9V_PDElMER4MwDaWyuCZ3KA5/view?usp=sharing），有什么办法可以加快优化的性能吗？我建模效率低下吗？如果我在代码中将“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. The 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()