-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathStarGenLite_SDUC.py
199 lines (166 loc) · 11.9 KB
/
StarGenLite_SDUC.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
# StarGen Lite Stochastic Daily Unit Commitment of Thermal and Hydro Units (SDUC)
# Developed by
# Andres Ramos
# Instituto de Investigacion Tecnologica
# Escuela Tecnica Superior de Ingenieria - ICAI
# UNIVERSIDAD PONTIFICIA COMILLAS
# Alberto Aguilera 23
# 28015 Madrid, Spain
# Andres.Ramos@comillas.edu
# MIT Energy Initiative
# Massachusetts Institute of Technology
# arght@mit.edu
# July 9, 2019
import numpy as np
import pandas as pd
from pyomo.environ import *
from pyomo.core import Var #To fix binary variables
from pyomo.core import Constraint #To get dual variables
from pyomo.core.kernel import *
# system dimensions
N = 24 # hours
SC = 3 # scenarios
G = 13 # thermal and hydro generating units
# reading data from Excel
InputFile = 'StarGenLite_SDUC.xlsm'
dfDemand = pd.read_excel(InputFile, sheet_name='DemandReserveIG', header=4, usecols='D' , nrows=N )
dfOperReserve = pd.read_excel(InputFile, sheet_name='DemandReserveIG', header=4, usecols='I' , nrows=N )
dfOperReserveUp = pd.read_excel(InputFile, sheet_name='DemandReserveIG', header=4, usecols='N' , nrows=N )
dfOperReserveDw = pd.read_excel(InputFile, sheet_name='DemandReserveIG', header=4, usecols='S' , nrows=N )
dfIntermGen = pd.read_excel(InputFile, sheet_name='DemandReserveIG', header=4, usecols='X:Z', nrows=N )
dfThermalGen = pd.read_excel(InputFile, sheet_name='Generation' , header=5, usecols='C:R', nrows=G, index_col=0)
# system parameters
pScenProb = [0.3, 0.5, 0.2] # probabilities of scenarios
pENSCost = 10 # cost of energy not served [MEUR per GWh]
pCO2Cost = 5 # cost of CO2 emission [EUR per tCO2]
#scaling of parameters to GW and MEUR
pDemand = dfDemand.values * 1e-3 # hourly load [GW]
pOperReserve = dfOperReserve.values * 1e-3 # hourly operating reserve [GW]
pOperReserveUp = dfOperReserveUp.values * 1e-3 # hourly operating reserve up [GW]
pOperReserveDw = dfOperReserveDw.values * 1e-3 # hourly operating reserve down [GW]
pIntermGen = dfIntermGen.values * 1e-3 # stochastic IG generation [GW]
pMaxProd = dfThermalGen.iloc[:, 0].values * 1e-3 # maximum output [GW]
pMinProd = dfThermalGen.iloc[:, 1].values * 1e-3 # minimum output [GW]
pIniOut = dfThermalGen.iloc[:, 2].values * 1e-3 # initial output [GW]
pRampUp = dfThermalGen.iloc[:, 3].values * 1e-3 # ramp up [GW/h]
pRampDw = dfThermalGen.iloc[:, 4].values * 1e-3 # ramp down [GW/h]
pSlopeVarCost = dfThermalGen.iloc[:, 6].values * 1e-3 * dfThermalGen.iloc[:, 5].values + dfThermalGen.iloc[:, 8].values * 1e-3 # slope variable cost [MEUR/GWh]
pInterVarCost = dfThermalGen.iloc[:, 7].values * 1e-6 * dfThermalGen.iloc[:, 5].values # intercept variable cost [MEUR/GWh]
pMinTU0 = dfThermalGen.iloc[:, 9].values # minimum time up [h]
pMinTD0 = dfThermalGen.iloc[:,10].values # minimum time down [h]
pEmissionCost = dfThermalGen.iloc[:,11].values * 1e-3 * pCO2Cost # emission cost [MEUR/GWh]
pStartUpCost = dfThermalGen.iloc[:,12].values * 1e-6 * dfThermalGen.iloc[:, 5].values # startup cost [MEUR]
pShutDownCost = dfThermalGen.iloc[:,13].values * 1e-6 * dfThermalGen.iloc[:, 5].values # shutdown cost [MEUR]
pIniUC = np.zeros(G)
pMinTU = np.zeros(G)
pMinTD = np.zeros(G)
for g in range(G):
if pIniOut[g] >= pMinProd[g]:
pIniUC[g] = 1
else:
pIniUC[g] = 0
if pMinTU0[g] == 0:
pMinTU[g] = 1
else:
pMinTU[g] = round(pMinTU0[g])
if pMinTD0[g] == 0:
pMinTD[g] = 1
else:
pMinTD[g] = round(pMinTD0[g])
# Stochastic Daily Unit Commitment (UC) SDUC
mSDUC = ConcreteModel()
# sets
mSDUC.N = Set(initialize=RangeSet(0,N -1), doc='hours' )
mSDUC.G = Set(initialize=RangeSet(0,G -1), doc='units' )
mSDUC.SC = Set(initialize=RangeSet(0,SC-1), doc='scenarios')
# parameters
pProduct1_up = np.zeros((SC,N,G))
pProduct_up = np.zeros((SC,N,G))
pIG_up = np.zeros((SC,N ))
pENS_up = np.zeros((SC,N ))
# definition of variable bounds
for sc in range(SC):
for n in range(N):
pIG_up [sc,n] = pIntermGen[n,sc]
pENS_up[sc,n] = pDemand[n]
for g in range(G):
pProduct1_up[sc,n,g] = pMaxProd[g]-pMinProd[g]
pProduct_up [sc,n,g] = pMaxProd[g]
# variable bounds
def vOutput_bd (mSDUC, sc,n,g):
return ((0.0, pProduct_up [sc,n,g]))
def vOutput2nd_bd(mSDUC, sc,n,g):
return ((0.0, pProduct1_up[sc,n,g]))
def vIG_bd (mSDUC, sc,n ):
return ((0.0, pIG_up [sc,n ]))
def vENS_bd (mSDUC, sc,n ):
return ((0.0, pENS_up [sc,n ]))
# variables
mSDUC.vOutput = Var(mSDUC.SC, mSDUC.N, mSDUC.G, within=NonNegativeReals, bounds=vOutput_bd , doc='output of the unit [GW]')
mSDUC.vOutput2nd = Var(mSDUC.SC, mSDUC.N, mSDUC.G, within=NonNegativeReals, bounds=vOutput2nd_bd, doc='output of the unit > min load [GW]')
mSDUC.vIG = Var(mSDUC.SC, mSDUC.N, within=NonNegativeReals, bounds=vIG_bd , doc='intermittent generation [GW]')
mSDUC.vENS = Var(mSDUC.SC, mSDUC.N, within=NonNegativeReals, bounds=vENS_bd , doc='energy not served [GW]')
mSDUC.vCommitment = Var( mSDUC.N, mSDUC.G, within=Binary , doc='commitment of the unit {0,1}')
mSDUC.vStartup = Var( mSDUC.N, mSDUC.G, within=Binary , doc='startup of the unit {0,1}')
mSDUC.vShutdown = Var( mSDUC.N, mSDUC.G, within=Binary , doc='shutdown of the unit {0,1}')
mSDUC.obj = Objective(expr=sum(pENSCost * mSDUC.vENS [sc,n ] * pScenProb[sc] for sc in mSDUC.SC for n in mSDUC.N ) +
sum(pSlopeVarCost[g] * mSDUC.vOutput [sc,n,g] * pScenProb[sc] for sc in mSDUC.SC for n in mSDUC.N for g in mSDUC.G) +
sum(pEmissionCost[g] * mSDUC.vOutput [sc,n,g] * pScenProb[sc] for sc in mSDUC.SC for n in mSDUC.N for g in mSDUC.G) +
sum(pInterVarCost[g] * mSDUC.vCommitment [ n,g] for n in mSDUC.N for g in mSDUC.G) +
sum(pStartUpCost [g] * mSDUC.vStartup [ n,g] for n in mSDUC.N for g in mSDUC.G) +
sum(pShutDownCost[g] * mSDUC.vShutdown[ n,g] for n in mSDUC.N for g in mSDUC.G) ,
sense=minimize, doc='total system variable cost [Meur]')
# constraints
def eBalance(mSDUC,sc,n):
return sum(mSDUC.vOutput[sc,n,g] for g in mSDUC.G) + mSDUC.vIG[sc,n] + mSDUC.vENS[sc,n] == float(pDemand[n])
mSDUC.eBalance = Constraint(mSDUC.SC, mSDUC.N, rule=eBalance, doc='load generation balance [GW]')
def eOpReserve(mSDUC,n):
return sum(pMaxProd[g] * mSDUC.vCommitment[n,g] for g in mSDUC.G) >= float(pOperReserve[n]) + float(pDemand[n])
mSDUC.eOpReserve = Constraint(mSDUC.N, rule=eOpReserve, doc='operating reserve [GW]')
def eReserveUp(mSDUC,sc,n):
return sum(pMaxProd[g] * mSDUC.vCommitment[n, g] - mSDUC.vOutput[sc,n,g] for g in mSDUC.G) >= float(pOperReserveUp[n])
mSDUC.eReserveUp = Constraint(mSDUC.SC, mSDUC.N, rule=eReserveUp, doc='operating reserve upwards [GW]')
def eReserveDw(mSDUC,sc,n):
return sum(pMinProd[g] * mSDUC.vCommitment[n,g] - mSDUC.vOutput[sc,n,g] for g in mSDUC.G) <= -float(pOperReserveDw[n])
mSDUC.eReserveDw = Constraint(mSDUC.SC, mSDUC.N, rule=eReserveDw, doc='operating reserve downwards [GW]')
def eMaxOutput(mSDUC,sc,n,g):
return mSDUC.vOutput[sc,n,g] / float(pMaxProd[g]) <= mSDUC.vCommitment[n,g]
mSDUC.eMaxOutput = Constraint(mSDUC.SC, mSDUC.N, mSDUC.G, rule=eMaxOutput, doc='max output of a committed unit [GW]')
def eMinOutput(mSDUC,sc,n,g):
return mSDUC.vOutput[sc,n,g] / float(pMinProd[g]) >= mSDUC.vCommitment[n,g]
mSDUC.eMinOutput = Constraint(mSDUC.SC, mSDUC.N, mSDUC.G, rule=eMinOutput, doc='min output of a committed unit [GW]')
def eTotOutput(mSDUC,sc,n,g):
return mSDUC.vOutput[sc,n,g] == float(pMinProd[g]) * mSDUC.vCommitment[n,g] + mSDUC.vOutput2nd[sc,n,g]
mSDUC.eTotOutput = Constraint(mSDUC.SC, mSDUC.N, mSDUC.G, rule=eTotOutput, doc='tot output of a committed unit [GW]')
def eRampUp(mSDUC,sc,n,g):
return mSDUC.vOutput2nd[sc,n,g] - max(float(pIniOut[g]-pMinProd[g]),0) <= float(pRampUp[g]) if n == 0 else mSDUC.vOutput2nd[sc,n,g] - mSDUC.vOutput2nd[sc,n-1,g] <= float(pRampUp[g])
mSDUC.eRampUp = Constraint(mSDUC.SC, mSDUC.N, mSDUC.G, rule=eRampUp, doc='bound on ramp up [GW]')
def eRampDw(mSDUC,sc,n,g):
return mSDUC.vOutput2nd[sc,n,g] - max(float(pIniOut[g]-pMinProd[g]),0) >= -float(pRampDw[g]) if n == 0 else mSDUC.vOutput2nd[sc,n,g] - mSDUC.vOutput2nd[sc,n-1,g] >= -float(pRampDw[g])
mSDUC.eRampDw = Constraint(mSDUC.SC, mSDUC.N, mSDUC.G, rule=eRampDw, doc='bound on ramp down [GW]')
def eUCStrShut(mSDUC,n,g):
return mSDUC.vCommitment[n,g] - float(pIniUC[g]) == mSDUC.vStartup[n,g] - mSDUC.vShutdown[n,g] if n == 0 else mSDUC.vCommitment[n,g] - mSDUC.vCommitment[n-1,g] == mSDUC.vStartup[n,g] - mSDUC.vShutdown[n,g]
mSDUC.eUCStrShut = Constraint(mSDUC.N, mSDUC.G, rule=eUCStrShut, doc='relation among commitment startup and shutdown')
def eMinTUp(mSDUC,n,g):
if n >= int(pMinTU[g]) & n <= N:
return sum(mSDUC.vStartup[nn,g] for nn in range(n+1-int(pMinTU[g]),n+1)) <= mSDUC.vCommitment[n, g]
else:
Constraint.Skip()
mSDUC.eMinTUp = Constraint(mSDUC.N, mSDUC.G, rule=eMinTUp, doc='minimum up time ( committed)')
def eMinTDw(mSDUC,n,g):
if n >= int(pMinTD[g]) & n <= N:
return sum(mSDUC.vShutdown[nn,g] for nn in range(n+1-int(pMinTD[g]),n+1)) <= 1 - mSDUC.vCommitment[n, g]
else:
Constraint.Skip()
mSDUC.eMinTDw = Constraint(mSDUC.N, mSDUC.G, rule=eMinTDw, doc='minimum down time (not committed)')
mSDUC.write('modelPyomo.lp', io_options={'symbolic_solver_labels': True}) # create file lp of the SDUC
solver = SolverFactory('gurobi') # Select solver
solver.options['MIPGap'] = 0.0 # Change Mipgap
results = solver.solve(mSDUC, tee=True) # tee=True displays the output of the solver
mSDUC.solutions.load_from(results) # necessary for fixing the values of the binary variables
# In case dual variables are required
mSDUC.dual = Suffix(direction=Suffix.IMPORT)
#pSRMC = np.zeros((SC,N))
#for sc in mSDUC.SC:
# for n in mSDUC.N:
# print('ppp', sc, n, mSDUC.dual[mSDUC.eBalance[sc,n]])