用于优化 EV 充电的非线性规划 APOPT 求解器在变量边界 <= 0 时并非不可行 (GEKKO python)
Nonlinear programming APOPT solver for optimal EV charging not infeasible with variables boundary <= 0 (GEKKO python)
我尝试使用 GEKKO 包进行最佳 EV 充电调度。但是，当我的代码设置为小于或等于零时，即 x=m.Array(m.Var,n_var,value=0,磅=0,ub=1.0）。错误信息是 'Unsuccessful with error code 0'。
下面是我的 python 脚本。如果您对此问题有任何建议，请随时告诉我。
谢谢，
奇猜
#------------------------------
import numpy as np
import pandas as pd
import math
import os
from gekko import GEKKO
print('...... Preparing data for optimal EV charging ........')
#---------- Read data from csv.file -----------)
Main_path = os.path.dirname(os.path.realpath(__file__))
Baseload_path = Main_path + '/baseload_data.csv'
TOU_path = Main_path + '/TOUprices.csv'
EV_path = Main_path + '/EVtravel.csv'
df_Baseload = pd.read_csv(Baseload_path, index_col= 'Time')
df_TOU = pd.read_csv(TOU_path, index_col= 'Time')
df_EVtravel = pd.read_csv(EV_path, index_col= 'EV_no')
#Create and change run directory
rd= r'.\RunDir'
if not os.path.isdir(os.path.abspath(rd)):
os.mkdir(os.path.abspath(rd))
#--------------------EV charging optimization function----------------#
def OptEV(tou_rate, EV_data, P_baseload):
"""This function is to run EV charging optimization for each houshold"""
#Initialize EV model and traveling data in 96 intervals(interval = 15min)
s_time= 4*(EV_data[0]//1) + math.ceil(100*(EV_data[0]%1)/15) #starting time
d_time= 4*(EV_data[1]//1) + math.floor(100*(EV_data[1]%1)/15) #departure time
Pch_rating= EV_data[2] #charing rated power(kW)
kWh_bat= EV_data[3] #Battery rated capacity(kWh)
int_SoC= EV_data[4] #Initial battery's SoC(p.u.)
#Calculation charging period
if d_time<= s_time:
ch_period = 96+d_time-s_time
else:
ch_period = d_time-s_time
Np= int(ch_period)
print('charging period = %d intervals'%(Np))
#Construct revelant data list based on charging period
ch_time = [0]*Np #charging time step list
price_rate = [0]*Np #electricity price list
kW_baseload = [0]*Np #kW house baseload power list
#Re-arrange charging time, electricity price rate and baseload
for i in range(Np):
t_step = int(s_time)+i
if t_step <= 95: #Before midnight
ch_time[i] = t_step
price_rate[i] = tou_rate[t_step]
kW_baseload[i] = P_baseload[t_step]/1000 #active house baseload
else: #After midnight
ch_time[i] = t_step-96
price_rate[i] = tou_rate[t_step-96]
kW_baseload[i] = P_baseload[t_step-96]/1000
#Initialize Model
m = GEKKO(remote=False) # or m = GEKKO() for solve locally
m.path = os.path.abspath(rd) # change run directory
#define parameter
ch_eff= m.Const(value=0.90) #charging/discharging efficiency
alpha= m.Const(value= 0.00005) #regularization constant battery profile
net_load= [None]*Np #net metering houshold load power array
elec_price= [None]*Np #purchased electricity price array
SoC= [None]*(Np+1) #SoC of batteries array
#initialize variables
n_var= Np #number of dicision variables
x = m.Array(m.Var,n_var,value=0,lb=0,ub=1.0) #dicision charging variables
#Calculation relevant evaluated parameters
#x[0] = m.Intermediate(-1.025)
SoC[0]= m.Intermediate(int_SoC) #initial battery SoC
for i in range(Np):
#Netload metering evaluation
net_load[i]= m.Intermediate(kW_baseload[i]+x[i]*Pch_rating)
#electricity cost price evaluation(cent/kWh)
Neg_pr= (1/4)*net_load[i]*price_rate[i] # Reverse power cost
Pos_pr= (1/4)*net_load[i]*price_rate[i] # Purchased power cost
elec_price[i]= m.Intermediate(m.if3(net_load[i], Neg_pr, Pos_pr))
#current battery's SoC evaluation
j=i+1
SoC_dch= (1/4)*(x[i]*Pch_rating/ch_eff)/kWh_bat #Discharging(V2G)
SoC_ch= (1/4)*ch_eff*x[i]*Pch_rating/kWh_bat #Discharging
SoC[j]= m.Intermediate(m.if3(x[i], SoC[j-1]+SoC_dch, SoC[j-1]+SoC_ch))
#m.solve(disp=False)
#-------Constraint functions--------#
#EV battery constraint
m.Equation(SoC[-1] >= 0.80) #required departure SoC (minimum=80%)
for i in range(Np):
j=i+1
m.Equation(SoC[j] >= 0.20) #lower SoC limit = 20%
for i in range(Np):
j=i+1
m.Equation(SoC[j] <= 0.95) #upper SoC limit = 95%
#household Net-power constraint
for i in range(Np):
m.Equation(net_load[i] >= -10.0) #Lower netload power limit
for i in range(Np):
m.Equation(net_load[i] <= 10.0) #Upper netload power limit
#Objective functions
elec_cost = m.Intermediate(m.sum(elec_price)) #electricity cost
#battery degradation cost
bat_cost = m.Intermediate(m.sum([alpha*xi**2 for xi in x]))
#bat_cost = 0 #Not consider battery degradation cost
m.Minimize(elec_cost + bat_cost) # objective
#Set global options
m.options.IMODE = 3 #steady state optimization
#Solve simulation
try:
m.solve(disp=True) # solve
print('--------------Results---------------')
print('Objective Function= ' + str(m.options.objfcnval))
print('x= ', x)
print('price_rate= ', price_rate)
print('net_load= ', net_load)
print('elec_price= ', elec_price)
print('SoC= ', SoC)
print('Charging time= ', ch_time)
except:
print('*******Not successful*******')
print('--------------No convergence---------------')
# from gekko.apm import get_file
# print(m._server)
# print(m._model_name)
# f = get_file(m._server,m._model_name,'infeasibilities.txt')
# f = f.decode().replace('\r','')
# with open('infeasibilities.txt', 'w') as fl:
# fl.write(str(f))
Pcharge = x
return ch_time, Pcharge
pass
#---------------------- Run scripts ---------------------#
TOU= df_TOU['Prices'] #electricity TOU prices rate (c/kWh)
Load1= df_Baseload['Load1']
EV_data = [17.15, 8.15, 3.3, 24, 0.50] #[start,final,kW_rate,kWh_bat,int_SoC]
OptEV(TOU, EV_data, Load1)
#--------------------- End of a script --------------------#
当求解器找不到解决方案并报告“找不到解决方案”时，有一个 troubleshooting method 来诊断问题。首先要做的是使用 m.solve(disp=True) 查看求解器输出。求解器可能已经确定了一个不可行的问题，或者它达到了最大迭代次数而没有收敛到一个解决方案。在你的情况下，它确定问题是不可行的。
不可行问题
如果求解器因方程不可行而失败，则它发现变量和方程的组合不可解。您可以尝试放宽变量边界或使用 运行 目录中的 infeasibilities.txt 文件确定哪个方程不可行。从本地 运行 目录中检索 infeasibilities.txt 文件，您可以在 m=GEKKO(remote=False).
时使用 m.open_folder() 查看该文件
最大迭代限制
如果求解器达到默认迭代限制 (m.options.MAX_ITER=250)，那么您可以尝试增加此限制或尝试以下策略。
- 通过为 APOPT 设置
m.options.SOLVER=1、为 BPOPT 设置 m.options.SOLVER=2、为 IPOPT 设置 m.options.SOLVER=3 或 m.options.SOLVER=0 来尝试不同的求解器来尝试所有可用的求解器。
- 首先求解变量数等于方程数的平方问题，求可行解。 Gekko 有几个选项可以帮助解决这个问题，包括
m.options.COLDSTART=1（为所有 FV 和 MV 设置 STATUS=0）或 m.options.COLDSTART=2（设置 STATUS=0 并执行块对角三角分解以找到可能的不可行方程）。
- 找到可行的解决方案后，尝试使用此解决方案作为初始猜测进行优化。
我尝试使用 GEKKO 包进行最佳 EV 充电调度。但是，当我的代码设置为小于或等于零时，即 x=m.Array(m.Var,n_var,value=0,磅=0,ub=1.0）。错误信息是 'Unsuccessful with error code 0'。 下面是我的 python 脚本。如果您对此问题有任何建议，请随时告诉我。
谢谢，
奇猜
#------------------------------
import numpy as np
import pandas as pd
import math
import os
from gekko import GEKKO
print('...... Preparing data for optimal EV charging ........')
#---------- Read data from csv.file -----------)
Main_path = os.path.dirname(os.path.realpath(__file__))
Baseload_path = Main_path + '/baseload_data.csv'
TOU_path = Main_path + '/TOUprices.csv'
EV_path = Main_path + '/EVtravel.csv'
df_Baseload = pd.read_csv(Baseload_path, index_col= 'Time')
df_TOU = pd.read_csv(TOU_path, index_col= 'Time')
df_EVtravel = pd.read_csv(EV_path, index_col= 'EV_no')
#Create and change run directory
rd= r'.\RunDir'
if not os.path.isdir(os.path.abspath(rd)):
os.mkdir(os.path.abspath(rd))
#--------------------EV charging optimization function----------------#
def OptEV(tou_rate, EV_data, P_baseload):
"""This function is to run EV charging optimization for each houshold"""
#Initialize EV model and traveling data in 96 intervals(interval = 15min)
s_time= 4*(EV_data[0]//1) + math.ceil(100*(EV_data[0]%1)/15) #starting time
d_time= 4*(EV_data[1]//1) + math.floor(100*(EV_data[1]%1)/15) #departure time
Pch_rating= EV_data[2] #charing rated power(kW)
kWh_bat= EV_data[3] #Battery rated capacity(kWh)
int_SoC= EV_data[4] #Initial battery's SoC(p.u.)
#Calculation charging period
if d_time<= s_time:
ch_period = 96+d_time-s_time
else:
ch_period = d_time-s_time
Np= int(ch_period)
print('charging period = %d intervals'%(Np))
#Construct revelant data list based on charging period
ch_time = [0]*Np #charging time step list
price_rate = [0]*Np #electricity price list
kW_baseload = [0]*Np #kW house baseload power list
#Re-arrange charging time, electricity price rate and baseload
for i in range(Np):
t_step = int(s_time)+i
if t_step <= 95: #Before midnight
ch_time[i] = t_step
price_rate[i] = tou_rate[t_step]
kW_baseload[i] = P_baseload[t_step]/1000 #active house baseload
else: #After midnight
ch_time[i] = t_step-96
price_rate[i] = tou_rate[t_step-96]
kW_baseload[i] = P_baseload[t_step-96]/1000
#Initialize Model
m = GEKKO(remote=False) # or m = GEKKO() for solve locally
m.path = os.path.abspath(rd) # change run directory
#define parameter
ch_eff= m.Const(value=0.90) #charging/discharging efficiency
alpha= m.Const(value= 0.00005) #regularization constant battery profile
net_load= [None]*Np #net metering houshold load power array
elec_price= [None]*Np #purchased electricity price array
SoC= [None]*(Np+1) #SoC of batteries array
#initialize variables
n_var= Np #number of dicision variables
x = m.Array(m.Var,n_var,value=0,lb=0,ub=1.0) #dicision charging variables
#Calculation relevant evaluated parameters
#x[0] = m.Intermediate(-1.025)
SoC[0]= m.Intermediate(int_SoC) #initial battery SoC
for i in range(Np):
#Netload metering evaluation
net_load[i]= m.Intermediate(kW_baseload[i]+x[i]*Pch_rating)
#electricity cost price evaluation(cent/kWh)
Neg_pr= (1/4)*net_load[i]*price_rate[i] # Reverse power cost
Pos_pr= (1/4)*net_load[i]*price_rate[i] # Purchased power cost
elec_price[i]= m.Intermediate(m.if3(net_load[i], Neg_pr, Pos_pr))
#current battery's SoC evaluation
j=i+1
SoC_dch= (1/4)*(x[i]*Pch_rating/ch_eff)/kWh_bat #Discharging(V2G)
SoC_ch= (1/4)*ch_eff*x[i]*Pch_rating/kWh_bat #Discharging
SoC[j]= m.Intermediate(m.if3(x[i], SoC[j-1]+SoC_dch, SoC[j-1]+SoC_ch))
#m.solve(disp=False)
#-------Constraint functions--------#
#EV battery constraint
m.Equation(SoC[-1] >= 0.80) #required departure SoC (minimum=80%)
for i in range(Np):
j=i+1
m.Equation(SoC[j] >= 0.20) #lower SoC limit = 20%
for i in range(Np):
j=i+1
m.Equation(SoC[j] <= 0.95) #upper SoC limit = 95%
#household Net-power constraint
for i in range(Np):
m.Equation(net_load[i] >= -10.0) #Lower netload power limit
for i in range(Np):
m.Equation(net_load[i] <= 10.0) #Upper netload power limit
#Objective functions
elec_cost = m.Intermediate(m.sum(elec_price)) #electricity cost
#battery degradation cost
bat_cost = m.Intermediate(m.sum([alpha*xi**2 for xi in x]))
#bat_cost = 0 #Not consider battery degradation cost
m.Minimize(elec_cost + bat_cost) # objective
#Set global options
m.options.IMODE = 3 #steady state optimization
#Solve simulation
try:
m.solve(disp=True) # solve
print('--------------Results---------------')
print('Objective Function= ' + str(m.options.objfcnval))
print('x= ', x)
print('price_rate= ', price_rate)
print('net_load= ', net_load)
print('elec_price= ', elec_price)
print('SoC= ', SoC)
print('Charging time= ', ch_time)
except:
print('*******Not successful*******')
print('--------------No convergence---------------')
# from gekko.apm import get_file
# print(m._server)
# print(m._model_name)
# f = get_file(m._server,m._model_name,'infeasibilities.txt')
# f = f.decode().replace('\r','')
# with open('infeasibilities.txt', 'w') as fl:
# fl.write(str(f))
Pcharge = x
return ch_time, Pcharge
pass
#---------------------- Run scripts ---------------------#
TOU= df_TOU['Prices'] #electricity TOU prices rate (c/kWh)
Load1= df_Baseload['Load1']
EV_data = [17.15, 8.15, 3.3, 24, 0.50] #[start,final,kW_rate,kWh_bat,int_SoC]
OptEV(TOU, EV_data, Load1)
#--------------------- End of a script --------------------#
当求解器找不到解决方案并报告“找不到解决方案”时，有一个 troubleshooting method 来诊断问题。首先要做的是使用 m.solve(disp=True) 查看求解器输出。求解器可能已经确定了一个不可行的问题，或者它达到了最大迭代次数而没有收敛到一个解决方案。在你的情况下，它确定问题是不可行的。
不可行问题
如果求解器因方程不可行而失败，则它发现变量和方程的组合不可解。您可以尝试放宽变量边界或使用 运行 目录中的 infeasibilities.txt 文件确定哪个方程不可行。从本地 运行 目录中检索 infeasibilities.txt 文件，您可以在 m=GEKKO(remote=False).
m.open_folder() 查看该文件
最大迭代限制
如果求解器达到默认迭代限制 (m.options.MAX_ITER=250)，那么您可以尝试增加此限制或尝试以下策略。
- 通过为 APOPT 设置
m.options.SOLVER=1、为 BPOPT 设置m.options.SOLVER=2、为 IPOPT 设置m.options.SOLVER=3或m.options.SOLVER=0来尝试不同的求解器来尝试所有可用的求解器。 - 首先求解变量数等于方程数的平方问题，求可行解。 Gekko 有几个选项可以帮助解决这个问题，包括
m.options.COLDSTART=1（为所有 FV 和 MV 设置 STATUS=0）或m.options.COLDSTART=2（设置 STATUS=0 并执行块对角三角分解以找到可能的不可行方程）。 - 找到可行的解决方案后，尝试使用此解决方案作为初始猜测进行优化。