Skip to content

Commit

Permalink
Merge pull request #43 from SutubraResearch/BugFixes
Browse files Browse the repository at this point in the history 
Implement a number of bug fixes
  • Loading branch information
@SutubraResearch
SutubraResearch authored Aug 18, 2023
2 parents 9063406 + 1b9442f commit b424dda
Show file tree
Hide file tree
Showing 4 changed files with 159 additions and 52 deletions.
55 changes: 28 additions & 27 deletions data_processing/DB_to_Excel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,10 +7,10 @@
from pyam import IamDataFrame

def make_excel(ifile, ofile, scenario):

if ifile is None :
raise "You did not specify the input file, remember to use '-i' option"
print("Use as :\n python DB_to_Excel.py -i <input_file> (Optional -o <output_excel_file_name_only>)\n Use -h for help.")
print("Use as :\n python DB_to_Excel.py -i <input_file> (Optional -o <output_excel_file_name_only>)\n Use -h for help.")
sys.exit(2)
else :
file_type = re.search(r"(\w+)\.(\w+)\b", ifile) # Extract the input filename and extension
Expand Down Expand Up @@ -75,18 +75,19 @@ def make_excel(ifile, ofile, scenario):
query = "SELECT regions, tech, sector, t_periods, emissions_comm, sum(emissions) as emissions FROM Output_Emissions WHERE scenario='" + scenario + "' GROUP BY \
regions, tech, sector, t_periods, emissions_comm"
df_emissions_raw = pd.read_sql_query(query, con)
df_emissions = df_emissions_raw.pivot_table(values='emissions', index=['regions', 'tech', 'sector','emissions_comm'], columns='t_periods')
df_emissions.reset_index(inplace=True)
df_emissions = pd.merge(all_emis_techs, df_emissions, on=['regions','tech', 'sector', 'emissions_comm'], how='left')
df_emissions.rename(columns={'regions':'Region', 'tech':'Technology', 'emissions_comm':'Emission Commodity', 'sector':'Sector'}, inplace=True)
df_emissions.to_excel(writer, sheet_name='Emissions', index=False, encoding='utf-8', startrow=1, header=False)
worksheet = writer.sheets['Emissions']
worksheet.set_column('A:A', 10)
worksheet.set_column('B:B', 10)
worksheet.set_column('C:C', 10)
worksheet.set_column('D:D', 20)
for col, val in enumerate(df_emissions.columns.values):
worksheet.write(0, col, val, header_format)
if not df_emissions_raw.empty:
df_emissions = df_emissions_raw.pivot_table(values='emissions', index=['regions', 'tech', 'sector','emissions_comm'], columns='t_periods')
df_emissions.reset_index(inplace=True)
df_emissions = pd.merge(all_emis_techs, df_emissions, on=['regions','tech', 'sector', 'emissions_comm'], how='left')
df_emissions.rename(columns={'regions':'Region', 'tech':'Technology', 'emissions_comm':'Emission Commodity', 'sector':'Sector'}, inplace=True)
df_emissions.to_excel(writer, sheet_name='Emissions', index=False, encoding='utf-8', startrow=1, header=False)
worksheet = writer.sheets['Emissions']
worksheet.set_column('A:A', 10)
worksheet.set_column('B:B', 10)
worksheet.set_column('C:C', 10)
worksheet.set_column('D:D', 20)
for col, val in enumerate(df_emissions.columns.values):
worksheet.write(0, col, val, header_format)

query = "SELECT regions, tech, sector, output_name, vintage, output_cost FROM Output_Costs WHERE output_name LIKE '%V_Discounted%' AND scenario='" + scenario + "'"
df_costs = pd.read_sql_query(query, con)
Expand Down Expand Up @@ -118,7 +119,7 @@ def make_excel(ifile, ofile, scenario):
df_activity['variable']='Activity|' + df_activity['sector'] + '|' + df_activity['tech']
df_activity.rename(columns={'t_periods':'year', 'vflow_out':'value', 'regions':'region'}, inplace=True)


# cast results to IamDataFrame and write to xlsx
columns = ['scenario', 'region', 'variable', 'year', 'value', 'unit']
_results = pd.concat([df_emissions_raw[columns], df_activity[columns], df_capacity[columns]])
Expand Down Expand Up @@ -147,10 +148,10 @@ def get_data(inputs):
ifile = None
ofile = None
scenario = set()

if inputs is None:
raise "no arguments found"

for opt, arg in inputs.items():
if opt in ("-i", "--input"):
ifile = arg
Expand All @@ -159,20 +160,20 @@ def get_data(inputs):
elif opt in ("-s", "--scenario"):
scenario.add(arg)
elif opt in ("-h", "--help") :
print("Use as :\n python DB_to_Excel.py -i <input_file> (Optional -o <output_excel_file_name_only>)\n Use -h for help.")
print("Use as :\n python DB_to_Excel.py -i <input_file> (Optional -o <output_excel_file_name_only>)\n Use -h for help.")
sys.exit()

make_excel(ifile, ofile, scenario)

if __name__ == "__main__":
if __name__ == "__main__":

try:
argv = sys.argv[1:]
opts, args = getopt.getopt(argv, "hi:o:s:", ["help", "input=", "output=", "scenario="])
except getopt.GetoptError:
print("Something's Wrong. Use as :\n python DB_to_Excel.py -i <input_file> (Optional -o <output_excel_file_name_only>)\n Use -h for help.")
sys.exit(2)
except getopt.GetoptError:
print("Something's Wrong. Use as :\n python DB_to_Excel.py -i <input_file> (Optional -o <output_excel_file_name_only>)\n Use -h for help.")
sys.exit(2)

print(opts)
get_data( dict(opts) )

get_data( dict(opts) )
2 changes: 1 addition & 1 deletion temoa_model/temoa_config.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,7 @@ def query_table (t_properties, f):
['param','SegFrac', '', '', 2],
['param','DemandSpecificDistribution','', '', 4],
['param','CapacityToActivity', '', '', 2],
['param','PlanningReserveMargin', '', '', 2],
['param','PlanningReserveMargin', '', '', 1],
['param','GlobalDiscountRate', '', '', 0],
['param','MyopicBaseyear', '', '', 0],
['param','DiscountRate', '', '', 3],
Expand Down
34 changes: 27 additions & 7 deletions temoa_model/temoa_initialize.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1016,14 +1016,34 @@ def DemandActivityConstraintIndices ( M ):
and "first_d" are the reference season and time-of-day, respectively used to
ensure demand activity remains consistent across time slices.
"""
first_s = M.time_season.first()
first_d = M.time_of_day.first()

# Find the first timestep of the year where the demand is non-zero:
eps = 0.000001
for r,p,t,v,dem in M.ProcessInputsByOutput.keys():
if dem in M.commodity_demand and t not in M.tech_annual:
for s in M.time_season:
for d in M.time_of_day:
if s != first_s or d != first_d:
yield (r,p,s,d,t,v,dem,first_s,first_d)
# No need for constraint if dem is not a demand commodity or
# if t is in tech_annual.
if (dem not in M.commodity_demand) or (t in M.tech_annual):
continue

# Find the first time step where the DSD is not 0
# Set the time_season and time_of_day to s0 and d0.
for s0 in M.time_season:
for d0 in M.time_of_day:
if (r,s0,d0,dem) in M.DemandSpecificDistribution.sparse_keys():
try:
if value(M.DemandSpecificDistribution[r,s0,d0,dem]) > eps:
break
except:
continue
else:
continue
break

# Start yielding the constraint indices
for s in M.time_season:
for d in M.time_of_day:
if s != s0 or d != d0:
yield (r,p,s,d,t,v,dem,s0,d0)

def DemandConstraintIndices ( M ):
used_dems = set((r,dem) for r, p, dem in M.Demand.sparse_iterkeys())
Expand Down
120 changes: 103 additions & 17 deletions temoa_model/temoa_rules.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -603,7 +603,7 @@ def CommodityBalance_Constraint(M, r, p, s, d, c):
interregional_exports = 0
if (r, p, c) in M.exportRegions:
interregional_exports = sum(
M.V_FlowOut[r+"-"+reg, p, s, d, c, S_t, S_v, S_o]
M.V_FlowOut[r+"-"+reg, p, s, d, c, S_t, S_v, S_o] / value(M.Efficiency[r+"-"+reg, c, S_t, S_v, S_o])
for reg, S_t, S_v, S_o in M.exportRegions[r, p, c]
)

Expand Down Expand Up @@ -696,7 +696,7 @@ def CommodityBalanceAnnual_Constraint(M, r, p, c):
interregional_exports = 0
if (r, p, c) in M.exportRegions:
interregional_exports = sum(
M.V_FlowOutAnnual[str(r)+"-"+str(reg), p, c, S_t, S_v, S_o]
M.V_FlowOutAnnual[str(r)+"-"+str(reg), p, c, S_t, S_v, S_o] / value(M.Efficiency[str(r)+"-"+str(reg), c, S_t, S_v, S_o])
for reg, S_t, S_v, S_o in M.exportRegions[r, p, c]
)

Expand Down Expand Up @@ -1446,30 +1446,48 @@ def ReserveMargin_Constraint(M, r, p, s, d):
r"""
During each period :math:`p`, the sum of the available capacity of all reserve
technologies :math:`\sum_{t \in T^{res}} \textbf{CAPAVL}_{r,p,t}`, which are
defined in the set :math:`\textbf{T}^{res}`, should exceed the peak load by
:math:`PRM`, the region-specific planning reserve margin. Note that the reserve
technologies :math:`\sum_{t \in T^{e}} \textbf{CAPAVL}_{r,p,t}`, which are
defined in the set :math:`\textbf{T}^{r,e}`, should exceed the peak load by
:math:`PRM`, the regional reserve margin. Note that the reserve
margin is expressed in percentage of the peak load. Generally speaking, in
a database we may not know the peak demand before running the model, therefore,
we write this equation for all the time-slices defined in the database in each region.
.. math::
:label: reserve_margin
\sum_{t \in T^{res}} {
\sum_{t \in T^{res} \setminus T^{e}} {
CC_{t,r} \cdot
\textbf{CAPAVL}_{p,t} \cdot
SEG_{s^*,d^*} \cdot C2A_{r,t} }
SEG_{s^*,d^*} \cdot C2A_{r,t} } +
\sum_{t \in T^{res} \cap T^{e}} {
CC_{t,r_i-r} \cdot
\textbf{CAPAVL}_{p,t} \cdot
SEG_{s^*,d^*} \cdot C2A_{r_i-r,t} } -
\sum_{t \in T^{res} \cap T^{e}} {
CC_{t,r-r_i} \cdot
\textbf{CAPAVL}_{p,t} \cdot
SEG_{s^*,d^*} \cdot C2A_{r_i-r,t} }
\geq
\sum_{ t \in T^{r,e},V,I,O } {
\textbf{FO}_{r, p, s, d, i, t, v, o} \cdot (1 + PRM_r)
{
\sum_{ t \in T^{res} \setminus T^{e},V,I,O }
\textbf{FO}_{r, p, s, d, i, t, v, o} +
\sum_{ t \in T^{res} \cap T^{e},V,I,O }
\textbf{FO}_{r_i-r, p, s, d, i, t, v, o} -
\sum_{ t \in T^{res} \cap T^{e},V,I,O }
\textbf{FI}_{r-r_i, p, s, d, i, t, v, o} -
\sum_{ t \in T^{res} \cap T^{s},V,I,O }
\textbf{FI}_{r, p, s, d, i, t, v, o}
}
\cdot (1 + PRM_r)
}
\\
\forall \{r, p, s, d\} \in \Theta_{\text{ReserveMargin}}
\text{and} \forall r_i \in R
"""
if (not M.tech_reserve) or ((r,p) not in M.processReservePeriods.keys()): # If reserve set empty or if r,p not in M.processReservePeriod.keys(), skip the constraint
if (not M.tech_reserve) or ((r, p) not in M.processReservePeriods.keys()): # If reserve set empty or if r,p not in M.processReservePeriod.keys(), skip the constraint
return Constraint.Skip

cap_avail = sum(
Expand All @@ -1482,22 +1500,90 @@ def ReserveMargin_Constraint(M, r, p, s, d):
if (r, p, t) in M.processVintages.keys()
for v in M.processVintages[r, p, t]
# Make sure (r,p,t,v) combinations are defined
if (r,p,t,v) in M.activeCapacityAvailable_rptv


)
if (r, p, t, v) in M.activeCapacityAvailable_rptv
)

# The above code does not consider exchange techs, e.g. electricity
# transmission between two distinct regions.
# We take exchange takes into account below.
# Note that a single exchange tech linking regions Ri and Rj is twice
# defined: once for region "Ri-Rj" and once for region "Rj-Ri".

# First, determine the amount of firm capacity each exchange tech
# contributes.
for r1r2 in M.RegionalIndices:
if '-' not in r1r2:
continue
r1, r2 = r1r2.split("-")

# Only consider the capacity of technologies that import to
# the region in question -- i.e. for cases where r2 == r.
if r2 != r:
continue

# add the available capacity of the exchange tech.
cap_avail += sum(
value(M.CapacityCredit[r1r2, p, t, v])
* M.ProcessLifeFrac[r1r2, p, t, v]
* M.V_Capacity[r1r2, t, v]
* value(M.CapacityToActivity[r1r2, t])
* value(M.SegFrac[s, d])
for t in M.tech_reserve
if (r1r2, p, t) in M.processVintages.keys()
for v in M.processVintages[r1r2, p, t]
# Make sure (r,p,t,v) combinations are defined
if (r1r2, p, t, v) in M.activeCapacityAvailable_rptv
)

# In most Temoa input databases, demand is endogenous, so we use electricity
# generation instead.
# generation instead as a proxy for electricity demand.
total_generation = sum(
M.V_FlowOut[r, p, s, d, S_i, t, S_v, S_o]
for (t,S_v) in M.processReservePeriods[r, p]
for (t, S_v) in M.processReservePeriods[r, p]
for S_i in M.processInputs[r, p, t, S_v]
for S_o in M.ProcessOutputsByInput[r, p, t, S_v, S_i]
)

cap_target = total_generation * (1 + value(M.PlanningReserveMargin[r]))
# We must take into account flows into storage technologies.
# Flows into storage technologies need to be subtracted from the
# load calculation.
total_generation -= sum(
M.V_FlowIn[r, p, s, d, S_i, t, S_v, S_o]
for (t, S_v) in M.processReservePeriods[r, p]
if t in M.tech_storage
for S_i in M.processInputs[r, p, t, S_v]
for S_o in M.ProcessOutputsByInput[r, p, t, S_v, S_i]
)

# Electricity imports and exports via exchange techs are accounted
# for below:
for r1r2 in M.RegionalIndices: # ensure the region is of the form r1-r2
if '-' not in r1r2:
continue
if (r1r2, p) not in M.processReservePeriods: # ensure the technology in question exists
continue
r1, r2 = r1r2.split("-")
# First, determine the exports, and subtract this value from the
# total generation.
if r1 == r:
total_generation -= sum(
M.V_FlowOut[r1r2, p, s, d, S_i, t, S_v, S_o] / value(M.Efficiency[r1r2, S_i, t, S_v, S_o])
for (t, S_v) in M.processReservePeriods[r1r2, p]
for S_i in M.processInputs[r1r2, p, t, S_v]
for S_o in M.ProcessOutputsByInput[r1r2, p, t, S_v, S_i]
)
# Second, determine the imports, and add this value from the
# total generation.
elif r2 == r:
total_generation += sum(
M.V_FlowOut[r1r2, p, s, d, S_i, t, S_v, S_o]
for (t, S_v) in M.processReservePeriods[r1r2, p]
for S_i in M.processInputs[r1r2, p, t, S_v]
for S_o in M.ProcessOutputsByInput[r1r2, p, t, S_v, S_i]
if (t, S_v) in M.processReservePeriods[r1r2, p]
)

cap_target = total_generation * (1 + value(M.PlanningReserveMargin[r]))
return cap_avail >= cap_target


Expand Down

0 comments on commit b424dda

Please sign in to comment.