# Simulation of SMDAMAGE, by John F Raffensperger. 10 Aug 2019, 10 Sep 2022, 30 Apr 2023, 17 May 2023. # All units should match Hector as defined in its input .ini files. # Input files, all from "Sources of data, corrected forestry.xlsm": # ../data/Forestry_Sequestration.csv, Year, Tonnes/hectare/year Loblolly pine, Tonnes/hectare/year Ponderosa pine Tonnes/hectare/year Black walnut # ../data/Forestry_bid_steps.csv, Tonnes/hectare/year Loblolly pine, Tonnes/hectare/year Ponderosa pine Tonnes/hectare/year Black walnut # ../data/MtC_bid_steps.csv, # ../data/C2F6_CF4_HFC125_HFC134a_HFC143a_SF6_bidsteps.csv, # ../data/CH4_bid_steps.csv, # ../data/N2O_bid_steps.csv. # ../data/Agriculture_bids.csv import matplotlib.pyplot as mplot from pulp import * #pulp.pulpTestAll() import time startTime = time.time() # 2020 <--------- bidding -----------2075 --------------------------------------------------------> 2120 # <--------- Temperature constrained ---------------------------------------------------------------------> 2301 # Timeline: StartYear, StartYear+1, ..., FirstConstrainedYear, ..., StartYear + 100, ..., FirstConstrainedYear + Number_of_constrained_years. # The whole timeline runs the PulseDataLength. StartYear = 2020 PeriodsPerYear = 1 PulseDataLength = 281 # Pulse data from Hector goes only 281 years. So raising this would understates later warming. LastConstraintYear = StartYear + PulseDataLength # Last period in the model, a constrained year. Number_of_bid_years = 150 BidPeriods = [float(StartYear) + float(t)/float(PeriodsPerYear) for t in range(PeriodsPerYear*Number_of_bid_years)] ModelPeriods = [float(StartYear) + float(t)/float(PeriodsPerYear) for t in range(PeriodsPerYear*PulseDataLength + 1)] # All objective function coefficients should be millions of dollars. So a bid of 1 is a bid for $1 million per unit of the chemical. # Caution, Carbon is 'ffi' in pulsefile. Chemicals = ['Agriculture', 'Black_walnut', 'C2F6', 'CF4', 'CH4', 'Carbon', 'HFC125', 'HFC134a', 'HFC143a', 'Loblolly_pine', 'N2O', 'Ponderosa_pine', 'Seaweed', 'SF6'] # Agriculture and "Carbon" are in megatons of carbon (not CO2). Hector uses gigatons of carbon, so we need to convert Hector's gigatons warming effects to SMDAMAGE megatons decision variables and back again to Hector gigatons for validation. Units = {'Agriculture':'mtC', 'Black_walnut_150':'mhectares', 'Black_walnut_10':'mhectares', 'Black_walnut_55':'mhectares', 'C2F6':'kt', 'CF4':'kt', 'CH4':'mt', 'Carbon':'mtC', 'HFC125':'kt', 'HFC134a':'kt', 'HFC143a':'kt', 'Loblolly_pine_150':'mhectares', 'Loblolly_pine_10':'mhectares', 'Loblolly_pine_24':'mhectares', 'N2O':'mt', 'Ponderosa_pine_150':'mhectares', 'Ponderosa_pine_10':'mhectares', 'Ponderosa_pine_103':'mhectares', 'Seaweed':'mt', 'SF6':'kt'} Treetypes = ['Loblolly_pine_150', 'Ponderosa_pine_150', 'Black_walnut_150', 'Loblolly_pine_10', 'Ponderosa_pine_10', 'Black_walnut_10', 'Loblolly_pine_24', 'Ponderosa_pine_103', 'Black_walnut_55', ] # 0. Preliminaries. -------------------------------------------------------------------------- Pulse = {} # [Pulseqty, warming1, warming2, warming3,...] # 1. Warming effects. ------------------------------------------------------------------------------------------ print ("1. Reading warming effects ", time.asctime(time.localtime(time.time()))) Wput_dict = {} # 1.1 Warming function. def Wput (pollutant, emissionperiod, warmingperiod): # Degrees Celsius/kg warming effect: pollutant p, emitted in period u, warming effect in period t. if (pollutant, emissionperiod, warmingperiod) in Wput_dict: return Wput_dict[(pollutant, emissionperiod, warmingperiod)] global Pulse global BidPeriods global Treetypes scaleCelsius = 1000.0 # Thousandths of a degree. assert (emissionperiod <= warmingperiod) assert(emissionperiod in BidPeriods) assert(warmingperiod in ModelPeriods) u_emissionperiod = BidPeriods.index(emissionperiod) t_warmingperiod = ModelPeriods.index(warmingperiod) if t_warmingperiod - u_emissionperiod >= PulseDataLength*PeriodsPerYear: return 0.0 # Units of Pulse input file are [ffi: GtC/yr, hence divide Carbon pulse by 1000, # CH4: MtCH4/yr, N2O: MtN2O-N/yr, C: MtC/yr, NMVOC: Mt/yr, BC: Mt/yr, OC: Mt/yr, # CF4: kt/yr, C2F6: kt/yr, HFC125: kt/yr, HFC134a: kt/yr, HFC143a: kt/yr, CFC11: kt/yr, CFC12: kt/yr, HCF22: kt/yr] # Agriculture. Degrees C in warmingperiod per million tons emitted in emissionperiod. Divide by 1000 because carbon pulse units are degrees C/gigaton. if pollutant == 'Agriculture': wput = -Pulse['ffi'][1 + t_warmingperiod - u_emissionperiod]*scaleCelsius/Pulse['ffi'][0]/1000.0 # Carbon. Degrees C in warmingperiod per million tons emitted in emissionperiod. Divide by 1000 because Carbon pulse units are degrees C/gigaton. elif pollutant == 'Carbon': # print (t_warmingperiod, u_emissionperiod, t_warmingperiod - u_emissionperiod) wput = Pulse['ffi'][1 + t_warmingperiod - u_emissionperiod]*scaleCelsius/Pulse['ffi'][0]/1000.0 # Seaweed. Degrees C in warmingperiod per million tons emitted in emissionperiod. Divide by 1000 because Carbon pulse units are degrees C/gigaton. elif pollutant == 'Seaweed': # print (t_warmingperiod, u_emissionperiod, t_warmingperiod - u_emissionperiod) wput = -Pulse['ffi'][1 + t_warmingperiod - u_emissionperiod]*scaleCelsius/Pulse['ffi'][0]/1000.0 # Forestry. Degrees C in warmingperiod per million tons sequestered in emissionperiod. Divide by 1000 because Carbon pulse units are degrees C/gigaton. elif pollutant in Treetypes: wput = Pulse[pollutant][1 + t_warmingperiod - u_emissionperiod]*scaleCelsius/Pulse[pollutant][0]/1000.0 if wput >= 0.0000001: # Forestry pulse should be negative. raise Exception("wput [" + pollutant + "," + str(u_emissionperiod) + "," + str(t_warmingperiod) + "] =" + str(wput) + " but should be negative.") else: # One of ['C2F6', 'CF4', 'CH4', 'HFC125', 'HFC134a', 'HFC143a', 'N2O', 'SF6'] if pollutant not in ['C2F6', 'CF4', 'CH4', 'HFC125', 'HFC134a', 'HFC143a', 'N2O', 'SF6']: raise Exception("Pollutant " + pollutant + " not in expected list.") wput = Pulse[pollutant][1 + t_warmingperiod - u_emissionperiod]*scaleCelsius/Pulse[pollutant][0] # Round off tiny values. Probably unnecessary. if abs(wput) < 10.0**-8: wput = 0.0 # if abs(wput) < 10.0**-6.5: wput = 0.0 Wput_dict[(pollutant, emissionperiod, warmingperiod)] = wput return wput # 1.2. Get warming effects for 'C2F6', 'CF4', 'CH4', 'Carbon', HFC125', 'HFC134a', 'HFC143a', 'N2O', 'SF6', 'SO2'. with open('../data/pulseoutput.txt', 'r') as pulsefile: # Annual warming. for line in pulsefile: # Each line looks like: ffi_emissions,13.931549999999998,GtC/yr,0.0,... chempulse = line.split(",") # Below, we're copying the annual warming for each period in the year. Pulse[chempulse[0].replace('_emissions','')] = [float(chempulse[1])] + [float(warming) for warming in chempulse[3:] for i in range(PeriodsPerYear)] # Units[chempulse[0]] = chempulse[2] # assert (Units[chempulse[0]] == chempulse[2]) # 1.3. Get carbon sequestered for forestry. Note that this is not warming. # For forestry, the year of planting is the "u_emissionperiod". Pulse units are kilograms, so we have to do the convolution. # We need the function Wput to calculate the Pulse value for Forestry. with open ('../data/Forestry_Sequestration.csv') as forestryfile: # Year, Tonnes/hectare/year Loblolly pine, Tonnes/hectare/year Ponderosa pine Tonnes/hectare/year Black walnut lines = [line.split(',') for line in forestryfile] # units = lines[0] tonnes = [] # convolution with carbon warming effect. for r, tree in enumerate(Treetypes): # Compare with the convolution in the spreadsheet. Pulse[tree] = [1.0] + [0.0 for t in ModelPeriods] for line in lines[1:]: # Year since planting, Tonnes/hectare/year Loblolly pine, Tonnes/hectare/year Ponderosa pine, Tonnes/hectare/year Black walnut # TonsSequesteredPerPeriod = float(line[r+1])/float(PeriodsPerYear) TonsSequesteredPerPeriod = float(line[r+1])/float(PeriodsPerYear) for p in range(PeriodsPerYear): sequesterperiod = float(line[0]) + float(p)/float(PeriodsPerYear) # Tree is planted in year 0. # Tree sequesters TonsSequesteredPerPeriod tonnes/hectare in each sequesterperiod from 0 to 155.75. # Decision variables are in millions of hectares/period, sequestration in millions of tons Carbon. # The pulse from sequesterperiod starts in sequesterperiod, and ends at ModelPeriods[-1]. # Now sum future periods of warming. sequesterperiodindex = ModelPeriods.index(StartYear + sequesterperiod) for ti, t in enumerate(ModelPeriods[sequesterperiodindex:]): Pulse[tree][1 + sequesterperiodindex + ti] -= float(TonsSequesteredPerPeriod)*Wput('Carbon', StartYear, t) # Examine warming. # Chemicals = ['Agriculture', 'Black_walnut', 'C2F6', 'CF4', 'CH4', 'Carbon', 'HFC125', 'HFC134a', 'HFC143a', 'Loblolly_pine', 'N2O', 'Ponderosa_pine', 'SF6'] # for c in Chemicals: # mplot.plot(ModelPeriods, [Wput(c,StartYear,t) for t in ModelPeriods]) # mplot.title(c) # mplot.show() # sys.exit() # 2. Bids. ------------------------------------------------------------------------------------------ print ("2. Reading bids...") # Parameters Bapt = {} # Bid price for agent a, pollutant p, time period t. Uapt = {} # Upper bid quantity, kg, for agent a, pollutant p, time period t. APT_set = set() # [(a, p, t),...] PT_set = set() # All bids in millions of dollars per million tons. # 2.1 Agriculture. $US2022, Mtons C2. Based on Smith P., et ak, 2014: Agriculture, Forestry and Other Land Use (AFOLU). In: Climate Change 2014: Mitigation of Climate Change, Figure 11.17. with open ('../data/Agriculture_bids.csv') as agfile: # Year, Tonnes/hectare/year Loblolly pine, Tonnes/hectare/year Ponderosa pine Tonnes/hectare/year Black walnut lines = [line.split(',') for line in agfile] Bid_Q_Agriculture = [(float (line[0]), float (line[1])) for line in lines[1:]] for t in BidPeriods: PT_set.add(('Agriculture',t)) for bidstep, bid in enumerate(Bid_Q_Agriculture): APT_set.add((bidstep,'Agriculture',t)) Bapt[bidstep,'Agriculture',t] = - bid[0] # M dollars/M tons Uapt[bidstep,'Agriculture',t] = bid[1]/float(PeriodsPerYear) # M tons per period. # 2.2 Carbon. $/tonne, Marginal Mtons Carbon. From "Sources of data.xlsm", sheet "Carbon bids", columns F,G. # Consider increasing initial quantity of 20 to 1000 or greater. with open('../data/MtC_bid_steps.csv', 'r') as Carbonfile: lines = [line for line in Carbonfile] for t in BidPeriods: PT_set.add(('Carbon',t)) for bidstep,line in enumerate(lines[1:]): bid = line.split(',') APT_set.add((bidstep,'Carbon',t)) Bapt[bidstep,'Carbon',t] = float(bid[0]) # m dollars/m tons Uapt[bidstep,'Carbon',t] = float(bid[1])/float(PeriodsPerYear) # Millions of tons bidstep += 1 # 2.3 Seaweed. $/tonne, Marginal Mtons Carbon. From "Sources of data.xlsm", sheet "Seaweed". with open('../data/Seaweed_bids.csv', 'r') as Seaweedfile: lines = [line for line in Seaweedfile] for t in BidPeriods: PT_set.add(('Seaweed',t)) for bidstep,line in enumerate(lines[1:]): bid = line.split(',') APT_set.add((bidstep,'Seaweed',t)) Bapt[bidstep,'Seaweed',t] = float(bid[0]) # m dollars/m tons Uapt[bidstep,'Seaweed',t] = float(bid[1])/float(PeriodsPerYear) # Millions of tons bidstep += 1 # Unused bid data: C6F14, HFC-152a, HFC-227ea, HFC-23, HFC245ca, HFC-32, HFC-43_10. # 2.3. Bid data retrieved here: C2F6, CF4, HFC-125, HFC-134a, HFC-143a, SF6. # File is sorted by chemical, year, price increasing. with open('../data/C2F6_CF4_HFC125_HFC134a_HFC143a_SF6_bidsteps.csv', 'r') as chemicalsfile: lines = [line for line in chemicalsfile] thislist = lines[0].split(',') # header: C2F6 $/kt,C2F6 kt/year,CF4 $/kt,CF4 kt/year,HFC125 $/kt,HFC125 kt/year,HFC134a $/kt,HFC134a kt/year,HFC143a $/kt,HFC143a kt/year,SF6 $/kt,SF6 kt/year chemicals = ['C2F6', 'CF4', 'HFC125', 'HFC134a', 'HFC143a', 'SF6'] # trust, but verify for chemical in chemicals: assert (chemical in Chemicals) for bidstep, line in enumerate(lines[1:]): # Skip header. thislist = line.split(',') thislist = [float(item) for item in thislist] for c, chemical in enumerate(chemicals): price, ktons = thislist[2*c:2*c+2] for t in BidPeriods: PT_set.add((chemicals[c],t)) APT_set.add((bidstep,chemicals[c],t)) # Data is given in $/tonne C. Quantities are kilotons, so decision variables should be $million/kiloton. # Thus, $500/ton --> $0.5 million per kiloton. Bapt[bidstep,chemicals[c],t] = round(price/1000.0, 6) # m dollars/k tons Uapt[bidstep,chemicals[c],t] = ktons/float(PeriodsPerYear) # Ktons per period. # 2.4. CH4_bid_steps.csv. with open('../data/CH4_bid_steps.csv', 'r') as chemicalsfile: chemical = 'CH4' lines = [line for line in chemicalsfile] for bidstep, line in enumerate(lines[1:]): price, mtons = line.split(',') price = float(price) mtons = float(mtons) for t in BidPeriods: PT_set.add((chemical,t)) APT_set.add((bidstep,chemical,t)) Bapt[bidstep,chemical,t] = price # Mega dollars/megatons # As described in "1-s2.0-S2352340919306882-mmc1, CH4 and NO2, jfr 1.xlsm", sheet "SSP2 CH4 N2O baseline emissions". Uapt[bidstep,chemical,t] = mtons/float(PeriodsPerYear) # 2.5. N2O_bid_steps.csv. with open('../data/N2O_bid_steps.csv', 'r') as chemicalsfile: chemical = 'N2O' lines = [line for line in chemicalsfile] for bidstep, line in enumerate(lines[1:]): price, mtons = line.split(',') price = float(price) mtons = float(mtons) for t in BidPeriods: PT_set.add((chemical,t)) APT_set.add((bidstep,chemical,t)) Bapt[bidstep,chemical,t] = price # Mega dollars/megatons # As described in "1-s2.0-S2352340919306882-mmc1, CH4 and NO2, jfr 1.xlsm", sheet "SSP2 CH4 N2O baseline emissions". Uapt[bidstep,chemical,t] = mtons/float(PeriodsPerYear) # Megatons # 2.6. Forestry_bid_steps.csv, with open ('../data/Forestry_bid_steps.csv') as forestryfile: # Tonnes/hectare/year Loblolly pine, Tonnes/hectare/year Ponderosa pine Tonnes/hectare/year Black walnut lines = [line.split(',') for line in forestryfile] # units = lines[0] # header: 'Plantable k hectares/year bid qty per crop', '2019 $/hectare Loblolly Pine', '2019 $/hectare ponderosa pine,2019 $/hectare black walnut'. bidstep = 0 for line in lines[1:]: # Plantable k hectares/year bid qty per crop, 2019 $/hectare Loblolly Pine, 2019 $/hectare ponderosa pine, 2019 $/hectare black walnut for t in BidPeriods: for r, tree in enumerate(Treetypes): PT_set.add((tree, t)) APT_set.add((bidstep, tree, t)) Bapt[bidstep, tree, t] = - float(line[1 + r]) # (M dollars)/(M hectares) Uapt[bidstep, tree, t] = float(line[0])/float(PeriodsPerYear) # M hectares plantable in each period. bidstep += 1 # ----------------------------------------------------------------------------------------------- Pollutants = sorted(list(set([p for p,t in PT_set]))) BidStepSet = {} for (p,t) in PT_set: BidStepSet[p,t] = [] TotalU = {(p,t): 0.0 for (p,t) in PT_set} for (a,p,t) in APT_set: TotalU[p,t] += Uapt[a,p,t] BidStepSet[p,t].append(a) # 3. Set up model. ------------------------------------------------------------------------------------------ ## Decision variables qapt = {(a,p,t): LpVariable("qapt(" + str(a) + "," + p + "," + str(t) + ")", 0.0, Uapt[a,p,t]) for (a, p, t) in APT_set} vpt = {(p,t): LpVariable("vpt(" + p + "," + str(t) + ")", None, None) for (p, t) in PT_set} # Must be a free variable. # for BeginConstraintYear in range(2075, 2076, 1): # range(2050, 2101, 1): for BeginConstraintYear in range(2110, 2180, 10): # range(2050, 2101, 1): # # Model print ("3. Creating model for " + str(BeginConstraintYear)) SMDAMAGE = LpProblem("SMDAMAGE", LpMaximize) # Objective ---------------------------------------------------------------------------------- print("Objective...") SMDAMAGE += lpSum(Bapt[a,p,t]*qapt[a,p,t] for (a, p, t) in APT_set), "Total value" # temperatureAboveSlack = {t: LpVariable("tempAboveSlack(" + str(t) + ")", 0, None) for t in ModelPeriods} # temperatureBelowSlack = {t: LpVariable("tempBelowSlack(" + str(t) + ")", 0, None) for t in ModelPeriods} # SMDAMAGE += lpSum ([-temperatureAboveSlack[t] + temperatureBelowSlack[t] for t in BidPeriods]), "Total temperature" # MAXIMIZING! # print (set([p for (a,p,t) in APT_set])) # should include all chemicals, forestry, agriculture. # # Constraint 3 in the paper. ------------------------------------------------------------- print("Vpt rows...") Vname = {} for (p, t) in PT_set: # SMDAMAGE += vpt[p,t] == lpSum (qapt[a,p1,t1] for (a,p1,t1) in APT_set if p1 == p and t1 == t), "Vpt(" + p + "," + str(t) + ")" Vname[(p,t)] = "Vpt(" + p + "," + str(t) + ")" SMDAMAGE += vpt[p,t] == lpSum ([qapt[a,p,t] for a in BidStepSet[p,t]]), Vname[(p,t)] # # Constraint 4 in the paper. ------------------------------------------------------------- print("Capt rows...") InitialBurden = {'C2F6': 62.0, 'CF4': 1125.1, 'CH4': 5153.0, 'Carbon': 248571.0, 'HFC125': 581.0, 'HFC134a': 1850.0, 'HFC143a': 352.0, 'N2O': 2546.0, 'SF6': 0.3} # Capt = {t: 200.0 if t < 2060 else 0.0 for t in ModelPeriods} # Allowed increase in global degrees Celsius in period t. # Capt = {t: 100.0 for t in ModelPeriods } # Allowed increase in global degrees Celsius in period t. # print (Capt) # exit() # Constrain each period. ConstraintPeriods = [float(BeginConstraintYear) + float(t)/float(PeriodsPerYear) for t in range(PeriodsPerYear*(PulseDataLength + StartYear - BeginConstraintYear))] Capt = {t: -1000.0 for t in ConstraintPeriods} # Allowed increase in global temp, thousandths of degrees Celsius in period t. # for t in ModelPeriods: SMDAMAGE += lpSum ([InitialBurden[p]*Wput(p,StartYear,t) for p in InitialBurden]\ # + [Wput(p,u,t)*vpt[p,u] for (p,u) in PT_set if u <= t])\ # + temperatureBelowSlack[t] - temperatureAboveSlack[t]\ # == Capt[t], "Capt(" + str(t) + ")" # for t in ModelPeriods: SMDAMAGE += lpSum ([InitialBurden[p]*Wput(p,StartYear,t) for p in InitialBurden]\ # + [Wput(p,u,t)*vpt[p,u] for (p,u) in PT_set if u <= t])\ # + temperatureBelowSlack[t] - temperatureAboveSlack[t]\ # == Capt[t], "Capt(" + str(t) + ")" temperatureChange = {t: LpVariable("tempChange(" + str(t) + ")", None, None) for t in ModelPeriods} for t in ModelPeriods: SMDAMAGE += lpSum ([InitialBurden[p]*Wput(p,StartYear,t) for p in InitialBurden]\ + [Wput(p,u,t)*vpt[p,u] for (p,u) in PT_set if u <= t]) - temperatureChange[t] == 0, "Temp_t(" + str(t) + ")" for t in ConstraintPeriods: SMDAMAGE += temperatureChange [t] <= Capt[t], "Capt(" + str(t) + ")" # Constrain every P periods. # temperatureSlack = {t: LpVariable("tempSlack(" + str(t) + ")", None, None) for t in ModelPeriods} # for t in ModelPeriods: SMDAMAGE += lpSum ([Wput(p,u,t)*vpt[p,u] for (p,u) in PT_set if u <= t]) - temperatureSlack[t] == Capt[t], "Capt(" + str(t) + ")" # print("Temp5 rows...") # sumperiod = 10 # ConstrainPyearPeriods = [t for t in ModelPeriods if t % sumperiod == 0] # Constraint5sumPeriods = {t:[t + float(u)/PeriodsPerYear for u in range(0, PeriodsPerYear*sumperiod)] for t in ConstrainPyearPeriods} # for t in ConstrainPyearPeriods: # SMDAMAGE += lpSum (temperatureSlack[u] for u in Constraint5sumPeriods[t]) <= 0.0, "Temp5Year(" + str(t) + ")" # print ("4. Writing a debug model...") ------------------------------------------------------------- # SMDAMAGE.writeLP(f"SMDAMAGE_{BeginConstraintYear}.lpt") print ("5. Solving the model for " + str(BeginConstraintYear)) solve_status = LpStatus[SMDAMAGE.solve(PULP_CBC_CMD(msg=0))] print (f"Done solving for begin constraint year {BeginConstraintYear}. Solve status: {solve_status}, Objective value = {value(SMDAMAGE.objective)}") # print ("6. Done solving. Solve status:", solve_status, "Objective value = ", value(SMDAMAGE.objective), time.asctime(time.localtime(time.time()))) # ReportPeriods = BidPeriods[0:4] + BidPeriods[100:104] + BidPeriods[-4:] # try: # Print to screen. # for p in Pollutants: # for t in ReportPeriods: print ("v[%s,%.2f] = %.1f" % (p, t, vpt[p, t].varValue)) # for p in Pollutants: # Show pi[p,t], dual prices for pollutants. # for t in ReportPeriods: print ("π[%s,%.2f] = $%.1f" % (p, t, SMDAMAGE.constraints[Vname[(p,t)]].pi)) # for t in ReportPeriods: # Show omega[t], dual prices for warming. # # print (t, SMDAMAGE.constraints["Capt(" + str(t) + ")"].pi, Capt[t]) # # if t in ModelPeriods[100:]: # print ("ω[%.2f] = $%.2f" % (t, SMDAMAGE.constraints["Capt(" + str(t) + ")"].pi)) # except: pass totalrevenue = 0.0 # Show net revenue with marginal cost pricing. for (p,t) in PT_set: totalrevenue -= vpt[p,t].varValue*SMDAMAGE.constraints[Vname[(p,t)]].pi print ("Total revenue = ", totalrevenue) # Save solution to CSV. with open ('smdamage_solution_' + str(BeginConstraintYear) + '.csv', 'w') as myoutputfile: myoutputfile.write("Solve status " + solve_status + ". Total revenue " + str(totalrevenue) + '\n') myoutputfile.write(','.join(['Year'] + [p + " " + Units[p] for p in Pollutants] + [p + " % max bid" for p in Pollutants] + [p + " $M/" + Units[p] for p in Pollutants]) + ',temp change' + ',Capt pi\n') for t in ModelPeriods[:-1]: line = [str(t)] # Year if t in BidPeriods: for p in Pollutants: line.append(str(vpt[p,t].varValue)) # qty units for p in Pollutants: line.append(str(vpt[p,t].varValue/TotalU[p,t])) # Fraction of bid for p in Pollutants: line.append(str(SMDAMAGE.constraints[Vname[(p,t)]].pi)) # $M/unit else: for p in Pollutants: line.append('\t\t\t') line.append(str(temperatureChange[t].varValue)) if t in ConstraintPeriods: line.append(str(SMDAMAGE.constraints["Capt(" + str(t) + ")"].pi)) else: line.append('.') myoutputfile.write (','.join(line) + '\n') # Tacks on multiple series within the for() loop. mplot.plot(ModelPeriods, [temperatureChange[t].varValue for t in ModelPeriods]) mplot.title(f"Change in global temp, 0.000C, 1st constraint year {BeginConstraintYear}") # mplot.show() myfigure = mplot.gcf() myfigure.savefig(f'Global_temp_first_constrained_in_{BeginConstraintYear}.svg', format='svg') myfigure.savefig(f'Global_temp_first_constrained_in_{BeginConstraintYear}.jpg', format='jpg') print ("Done, %s, %.1f seconds." % (time.asctime(time.localtime(time.time())), float(time.time() - startTime))) print ("Reminder: convert $/ton C to $/ton CO2.") # Graphs # Prices, all pollutants on one graph. # Log quantities, all pollutants on one graph. # Total warming by year.