I have tried to do optimal EV charging scheduling using the GEKKO packages. However, my code is stuck on some variable boundary condition when it is set to be lower than or equal to zero, i.e., x=m.Array(m.Var,n_var,value=0,lb=0,ub=1.0). The error message is 'Unsuccessful with error code 0'. Below is my python script. If you have any advice on this problem, please don't hesitate to let me know.
Thanks,
Chitchai
#------------------------------
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 --------------------#
When the solver fails to find a solution and reports "Solution Not Found", there is a troubleshooting method to diagnose the problem. The first thing to do is to look at the solver output with m.solve(disp=True). The solver may have identified either an infeasible problem or it reached the maximum number of iterations without converging to a solution. In your case, it identified the problem as infeasible.
Infeasible Problem
If the solver failed because of infeasible equations then it found that the combination of variables and equations is not solvable. You can try to relax the variable bounds or identify which equation is infeasible with the infeasibilities.txt file in the run directory. Retrieve the infeasibilities.txt file from the local run directory that you can view with m.open_folder() when m=GEKKO(remote=False).
Maximum Iteration Limit
If the solver reached the default iteration limit (m.options.MAX_ITER=250) then you can either try to increase this limit or else try the strategies below.
m.options.SOLVER=1 for APOPT, m.options.SOLVER=2 for BPOPT, m.options.SOLVER=3 for IPOPT, or m.options.SOLVER=0 to try all the available solvers.m.options.COLDSTART=1 (sets STATUS=0 for all FVs and MVs) or m.options.COLDSTART=2 (sets STATUS=0 and performs block diagonal triangular decomposition to find possible infeasible equations).