Commit ecb66c21 authored by ruzhu's avatar ruzhu
Browse files

Upload New File

parent 1c553e5f
Pipeline #21314 failed with stages
in 28 seconds
# -*- coding: utf-8 -*-
"""
Created on Mon Mar 15 11:42:56 2021
@author: ruzhu
"""
import pandas as pd
import numpy as np
from gams import *
import os
import sys
import subprocess
import rainflow
import math
import matplotlib.pyplot as plt
'''#Getting Started
# Step1: let python find GAMS
sys.path.append(r'C:\path\to\gams\apifiles\Python\api_37')
sys.path.append(r'C:\path\to\gams\apifiles\Python\gams')
# Step2: Copy the GAMS Python files to the Python installation
#Open a command prompt.
#cd C:\path\to\gams\apifiles\Python\api_37
#python setup.py install
#pip install rainflow
'''
def ReadData():
wind_forecasts=pd.DataFrame()
market_prices=pd.DataFrame()
for i in range(2016,2021):
market_price=pd.read_excel('./Data/elspot-prices_'+str(i)+'_hourly_eur.xlsx', engine='openpyxl')
wind_forecast=pd.read_excel('./Data/wind-power-dk-prognosis_'+str(i)+'_hourly.xlsx', engine='openpyxl')
wind_forecasts = pd.concat([wind_forecasts,wind_forecast['DK2']])
market_prices = pd.concat([market_prices,market_price['DK2']])
return wind_forecasts, market_prices
def Linear_Degfun(
kdelta1, kdelta2, kdelta3,
ksigma, sigma_ref, kT, Tref, kti,
rf_DoD, rf_SoC, rf_count, i
):
ld = 0
for j in range(0, len(rf_DoD)):
S_DoD = (kdelta1*rf_DoD[j]**kdelta2+kdelta3)**(-1)
S_time = kti*((i+1)*24/sum(rf_count)*rf_count[j]*3600)
S_SoC = math.exp(ksigma*(rf_SoC[j]-sigma_ref))
S_T = math.exp(kT*(20-Tref)*Tref/20)
ld = ld+(S_DoD+S_time)*S_SoC*S_T
return ld
def Non_Linear_Degfun(alpha, beta, ld, pre_nld, ld0):
if (pre_nld<0.8):
nld = 1-alpha*math.exp(-ld*beta)-(1-alpha)*math.exp(-ld)
else:
nld1 = 1-alpha*math.exp(-ld0*beta)-(1-alpha)*math.exp(-ld0)
nld = 1-(1-nld1)*math.exp(-(ld-ld0))
return nld
def Degcalculation(whole_SoC, ld0):
alpha = 0.0575
beta = 121
kdelta1 = 140000
kdelta2 = -0.5010
kdelta3 = -123000
ksigma = 1.04
sigma_ref = 0.5
kT = 0.0693
Tref = 25
kti = 4.14e-10
rf_DoD = list()
rf_SoC = list()
rf_count = list()
for rng, mean, count, i_start, i_end in rainflow.extract_cycles(whole_SoC):
rf_DoD.append(rng)
rf_SoC.append(mean)
rf_count.append(count)
ld = Linear_Degfun(
kdelta1, kdelta2, kdelta3,
ksigma, sigma_ref, kT, Tref, kti,
rf_DoD, rf_SoC, rf_count, i
)
nld = Non_Linear_Degfun(alpha, beta, ld, pre_nld, ld0)
return ld, nld, rf_DoD, rf_SoC, rf_count
def slope_update(whole_Pdis, whole_Pcha, whole_nld, pre_nld, nld, i):
if (i<8):
ad = (nld-pre_nld)/sum(whole_Pdis[-24:]+Pcha[-24:])
else:
ad = (whole_nld[-1]-whole_nld[-8])/sum(whole_Pdis[-24*7:]+Pcha[-24*7:])
return ad
if __name__ == "__main__":
#Call GAMS
if len(sys.argv) > 1:
ws = GamsWorkspace(system_directory = sys.argv[1])
else:
ws = GamsWorkspace(working_directory='GAMS')
#Initial data
ad = 0
nld = 0
ld0 = 0
pre_nld = 0
EBESS = 150
PwMax = 100
whole_SoC = list()
whole_Pdis = list()
whole_Pcha = list()
whole_nld = list()
whole_Profit = list()
whole_Deg_cost = list()
whole_Pw_cur = list()
wind_forecasts, market_prices = ReadData()
wind_forecasts = wind_forecasts/wind_forecasts.max()*PwMax
# Write GDX file for .gms
''' a = db.add_parameter("EBESS", 0)
a.add_record().value = EBESS
a = db.add_parameter("PwMax", 0)
a.add_record().value = PwMax'''
# for i in range(0,int(len(wind_forecasts)/24-1)):
for i in range(0,10):
# Write GDX file for .gms
# if os.path.exists(".\GAMS\DEMS_inputs.gdx"):
# os.remove(".\GAMS\DEMS_inputs.gdx")
db = ws.add_database()
iiset = list(range(1,25))
iiset = [str(x) for x in iiset]
ii = db.add_set("ii", 1)
for jj in iiset:
ii.add_record(jj)
wind_forecast=wind_forecasts[24*i:24*(i+1)]
wind_forecast=wind_forecast.values.tolist()
wind_forecast=[jj for ii in wind_forecast for jj in ii]
a = db.add_parameter_dc("Pw_hat", [ii])
for j in iiset:
a.add_record(j).value = wind_forecast[int(j)-1]
market_price=market_prices[24*i:24*(i+1)]
market_price=market_price.values.tolist()
market_price=[jj for ii in market_price for jj in ii]
a = db.add_parameter_dc("Sptprc_hat", [ii])
for j in iiset:
a.add_record(j).value = market_price[int(j)-1]
a = db.add_parameter("Emax", 0)
a.add_record().value = EBESS*(1-nld)
a = db.add_parameter("ad", 0)
a.add_record().value = ad
db.export("DEMS_inputs.gdx")
#Operate DEMS.gms
t1 = ws.add_job_from_file("DEMS.gms")
# opt = ws.add_options()
# opt.defines['suppress'] = "SuppressCompilerListing"
# opt.suppress = Suppress.SuppressCompilerListing
# opt.Listing = Listing.OffListing
t1.run(create_out_db = False)
# Read GDX file from .gms
db2 = ws.add_database_from_gdx("DEMS_outputs.gdx")
Pdis = {rec.keys[0]:rec.level for rec in db2["Pdis"]}
Pcha = {rec.keys[0]:rec.level for rec in db2["Pcha"]}
SoC = {rec.keys[0]:rec.level for rec in db2["SoC"]}
Profit = db2["obj"].first_record().level
Deg_cost = db2["Dcost"].first_record().level
Pw = {rec.keys[0]:rec.level for rec in db2["Pw"]}
Pdis = list(Pdis.values())
Pcha = list(Pcha.values())
SoC = list(SoC.values())
Pw = list(Pw.values())
whole_SoC = whole_SoC+SoC[0:24]
whole_Pdis = whole_Pdis+Pdis[0:24]
whole_Pcha = whole_Pcha+Pcha[0:24]
whole_Profit.append(Profit)
whole_Deg_cost.append(Deg_cost)
whole_Pw_cur = whole_Pw_cur+list(np.array(wind_forecast)-np.array(Pw[0:24]))
# Degradation Calculation
ld, nld, rf_DoD, rf_SoC, rf_count = Degcalculation(whole_SoC, ld0)
ad = slope_update(whole_Pdis, whole_Pcha, whole_nld, pre_nld, nld, i)
ld0 = ld
pre_nld = nld
whole_nld.append(nld)
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment