变量边界 <= 0 (GEKKO python) 的非线性规划 APOPT 求解器可实现最佳 EV 充电

Question

我已尝试使用 GEKKO 包进行最佳 EV 充电调度。但是，当我的代码设置为小于或等于零时，我的代码会卡在一些变量边界条件上，即 x=m.Array(m.Var,n_var,value=0,lb=0,ub=1.0) .错误消息是“不成功，错误代码 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 --------------------#

Answer 1

当求解器找不到解决方案并报告“未找到解决方案”时，有一个troubleshooting method 来诊断问题。首先要做的是使用m.solve(disp=True) 查看求解器输出。求解器可能已经确定了一个不可行的问题，或者它达到了最大迭代次数而没有收敛到一个解决方案。在您的情况下，它认为问题是不可行的。

不可行的问题

如果求解器由于不可行的方程而失败，那么它会发现变量和方程的组合是不可解的。您可以尝试放宽变量界限或使用运行目录中的infeasibilities.txt 文件确定哪个方程不可行。从本地运行目录中检索infeasibilities.txt 文件，当m=GEKKO(remote=False) 时，您可以使用m.open_folder() 查看该文件。

最大迭代限制

如果求解器达到默认迭代限制 (m.options.MAX_ITER=250)，那么您可以尝试增加此限制或尝试以下策略。

• 通过设置m.options.SOLVER=1 为APOPT、m.options.SOLVER=2 为BPOPT、m.options.SOLVER=3 为IPOPT 或m.options.SOLVER=0 来尝试不同的求解器以尝试所有可用的求解器。
• 首先通过求解变量数等于方程数的平方问题找到可行解。 Gekko 有几个选项可以帮助解决这个问题，包括 m.options.COLDSTART=1（为所有 FV 和 MV 设置 STATUS=0）或m.options.COLDSTART=2（设置 STATUS=0 并执行块对角三角分解以找到可能的不可行方程）。
• 一旦找到可行的解决方案，请尝试将此解决方案作为初始猜测进行优化。