%=====================================================
% Filename: ConstantsEM.m
% Copyright (c) 2015 John George Anasis
% **FULL THESIS MODEL WITH ANNUAL EMISSIONS AS THE PRIMARY CONSTRAINT.
% Linearly decreasing emissions limit to a final limit of 5.2Gt in 2025.
% No geoengineering other than tree planting assumed.
% **CASE WITH 50 MILLION HECTARES OF TREES PLANTED AT THE START OF THE
% SIMULATION TO SEE HOW THAT AFFECTS EMISSIONS.
% **NON-ENERGY CH4 AND N2O EMISSIONS REDUCED BY 15% OVER THE 10 YEAR
% SIMULATION.
% **CASE WITH PER UNIT OIL, GAS, AND COAL CO2 EMISSIONS UPDATED BASED ON
% EIA DATA.
% **REVISED LOWER BOUND COMPUTATION TO SET LOWER BOUND BASED ON AN
% ALLOWABLE DECREASE FROM THE RESOURCE'S VALUE AT THE PREVIOUS TIME STEP.
% **EXTENDED VERSION WITH UPDATED CO2 MODEL AND DETAILED ENERGY DEMANDS AND
% **VARIABLES TO ACCOUNT FOR SPECIFIC TRANSPORTATION AND INDUSTRIAL ENERGY
% **FOR COMBUSTABLE FUELS.
% **CASE WITH ENERGY EFFICIENCY CODE REVISED TO ALLOW FOR VARIABLE AMOUNT
% OF EFFICIENCY TO BE SPECIFIED.
% **CODE UPDATED SO THAT COSTS CAN VARY OVER TIME. CAN BE IMPLEMENTED
% EITHER AS A FUNCTION OR BY READING IN AN EXCEL SPREADSHEET. SEPARATED
% VARIABLE COSTS FOR COMBUSTION AND ELECTRICITY PRODUCTION FOR COAL AND
% GAS JUST AS WAS DONE FOR CAPITAL COSTS.
% **ADDED SEA SPRAY CLOUD ALBEDO ENHANCEMENT TO GEOENGINEERING OPTIONS.
% **CORRECTED SEA SPRAY EQUIVALENT EMISSIONS.
% **ASSUMED ZERO INITIAL FOREST AREA.
% **UPDATED CODE TO MAKE RESOURCE CHANGES MORE REALISTIC (GRADUAL PHASE OUT
% OF INITIAL RESOURCES, BETTER MANAGEMENT OF RESOURCE CAPS
% INCLUDING PHASE OUT OF OIL, GAS, AND COAL AS THESE RESOURCES ARE USED UP
% AND REALISTIC RESOURCE RAMP UPS).
%
% **ENERGY EFFICIENCY APPLIED TO LIQUID DEMAND.
% **COST OF ENERGY EFFICIENCY FIXED AS A PERCENTAGE OF STARTING YEAR GDP.
%
% **TREE PLANTING GEOENGINEERING EFFECTIVENESS BASED ON LOGISTIC (SIGMOID)
% FUNCTION RATHER THAN A SIMPLE EXPONENTIAL DECAY.
% **Added non-energy methane and N2O emissions in order to obtain more
% complete radiative forcing.
% **MODEL WITH CARBON TAX INCLUDED.
% **LOOKING AT THE WORLD AS A WHOLE (NO REGIONAL VARIATIONS IN PRICE OR
% **TECHNOLOGY).
%
% **RADIATIVE FORCING COMPUTED BASED ON PRE-INDUSTRIAL CO2 CONCENTRATION
% **RATHER THAN FROM THE STARTING YEAR CONCENTRATION.
% **COMPLETELY REVISED CAPITAL COST EQUATIONS. CAPITAL COSTS NOW MODELED AS
% AN ANNUAL VALUE (i.e. ANNUAL CAPITAL RECOVERY VALUE) AND COMBINED WITH
% THE VARIABLE & PENALTY COSTS IN THE OBJECTIVE FUNCTION AS ONE COST VALUE. ANY
% DECOMMISSIONING COSTS ALSO TURNED INTO ANNUAL COSTS AND COMBINED WITH THE
% ANNUAL CAPITAL COSTS.
% **Initial starting values are the results of the previous time step.
% **Initial x0 value at first time step set to resource values at t=0.
% **USE MATLAB OPTIMIZATION TOOLBOX fmincon FUNCTION TO DO THE OPTIMIZATION.
% REDUCED THE REQUIRED RESERVES TO 1% INSTEAD OF 10%
% REMOVED LOW TEMPARTURE CHANGE CONSTRAINT. TEMP CHANGE CAN NOW GO
% NEGATIVE.
% INCREASED VALUE OF MAXFUNEVALS PARAMETER TO 15,000 TO GET INTERATION 1 TO SOLVE.
% **MODEL INCLUDES THE FOLLOWING RESOURCE AND GEOENGINEERING OPTIONS
% WITH THE CORRESPONDING INDEX NUMBERS:
% 1 - oil(transport), 2 - natural gas(combustion), 3 - coal(combustion),
% 4 - nuclear, 5 - hydro, 6 - wind, 7 - geothermal, 8 - solar PV,
% 9 - solar thermal, 10 - biofuel(transport), 11 - biomass,
% 12 - energy efficiency, 13 - sulfur injection, 14 - iron seeding,
% 15 - tree planting, 16 - oil(reserve), 17 - natural gas(reserve),
% 18 - coal(reserve), 19 - hydro(reserve), 20 - geothermal(reserve),
% 21 - biofuel(reserve), 22 - biomass, 23 - sea spray,
% 24 - oil(combustion), 25 - gas(electric), 26 - coal(electric)
% 27 - coal(transport), 28 - biofuel(combustion)
% July 31, 2015
%=======================================================
function ConstantsEM
% All constants will be in a .mat file called "Plant".
% input length of time horizon (time), number of resources and geoengineering
% options.
time = 100;
Plant.time = time;
TotalResources = 28; % Need to include the extra variables needed for the reserves
Plant.TotalResources = TotalResources;
Resources = 21; % Total actual resources, not including any reserves.
Plant.Resources = Resources;
ReserveNumber = 7; % Total number of resources used for reserves.
Plant.ReserveNumber = ReserveNumber;
% define constants
Plant.TempLim = 2.0; % Maximum global mean temperature rise by end of the time horizon of study.
Plant.dis_rate = .04; % discount rate i (% shown as a decimal)
dis_rate = Plant.dis_rate; % Needed to convert capital costs to annual values.
Plant.tau1 = 394.4; % First lifetime parameter of CO2 (years)
Plant.tau2 = 36.54; % Second lifetime parameter of CO2 (years)
Plant.tau3 = 4.304; % Third lifetime parameter of CO2 (years)
Plant.a1 = 0.2173; % First coefficient for CO2 equation (unitless)
Plant.a2 = 0.2240; % Second coefficient for CO2 equation (unitless)
Plant.a3 = 0.2824; % Third coefficient for CO2 equation (unitless)
Plant.a4 = 0.2763; % Fourth coefficient for CO2 equation (unitless)
Plant.taum = 12.4; % lifetime of methane in the atmosphere (years)
Plant.taun = 121; % lifetime of N2O in the atmosphere (years)
Plant.GWPm = 34; % 100 year Global Warming Potential of methane.
Plant.GWPn = 268; % 100 year Global Warming Potential of N2O
Plant.GWPs = -58; % 100 year Global Warming Potential of SO2 injection.
Plant.co2i = 375; % initial atmospheric CO2 concentration at time t = 0 in ppm
Plant.ch4i = 1750; % initial atmospheric methane concentration at time t = 0 in ppb
Plant.n2oi = 320; % initial atmospheric N2O concentration at time t = 0 in ppb.
%Plant.n2oi = 270; % pre-industrial atmospheric N2O concentration in ppb.
Plant.R = 1E-4; % desired reserve margin over and above actual energy demand.
Plant.Lp = 3E3; % total acreage available for biofuel and biomass production plus tree planting (million ha).
% Specify the carbon tax. Convert tax in $/metric ton CO2-eq to
% $ trillion/gram CO2-eq emitted in order to line up with the other
% units used in the code.
CarbonTax = 0; % Carbon tax in $/metric ton CO2-eq.
Plant.ctax = CarbonTax * 1E-18; % Conversion to $trillion/gram CO2-eq.
% define the desired emissions limit.
% 2005 US emissions were 7.2Gt, so a 28% reduction by 2025 is a target
% of 5.18Gt. Year 2016 emissions from no limits scenario were
% 6.92Gt, so a reduction of 0.19Gt per year is necessary to reach the
% final target.
emtarget = 60; % Annual emissions target in billion metric tonnes CO2-eq (CDE) per year.
%emstep = -0.19; % Annual change in emissions target in billion metric tonnes CO2-eq.
%empercent = 0.98; % Annual change in emissions target as a fraction of the prior year target
Plant.emlimit = emtarget .* ones(time,1); % Matrix of emissions targets for constraint equation.
% Emissions targets matrix assuming annual targets are not fixed. Can
% choose two versions: A fixed change at each time step or a percentage
% at each time step.
% emlimit = zeros(time,1);
% for v = 1:time;
% if v == 1;
% emlimit(v) = emtarget;
% else
% emlimit(v) = emlimit(v-1) + emstep;
% %emlimit(v) = emlimit(v-1) * empercent;
% end
% end
% Plant.emlimit = emlimit;
%
% Oil Parameters (billions of barrels)
Plant.Co = 9.8E-2*ones(time,1); % variable cost of a barrel of oil ($trillion/ billion barrels)
%Plant.Co = xlsread('oilvarcost.xls'); % variable cost of a barrel of oil ($trillion/ billion barrels)
Co = Plant.Co;
%Plant.Ko = 0*ones(time,1); % capital cost of a barrel of oil ($trillion/ billion barrels)
Plant.Ko = 3.38E-2*ones(time,1); % capital cost of a barrel of oil ($trillion/ billion barrels)
%Plant.Ko = xlsread('oilcapcost.xls'); % capital cost of a barrel of oil ($trillion/ billion barrels)
Ko = Plant.Ko;
po = 0*ones(time,1); % penalty cost of a barrel of oil ($trillion/ billion barrels)
Plant.ko = 5.80; % conversion of billion bbl of oil to EJ (EJ/ billion bbl)
Plant.etao = 0.28; % efficiency of of burning oil for transportation needs.
Plant.etaoc = 0.72; % efficiency of of burning oil for specific industrial combustion/furnace needs.
Plant.Lo = 3.35E3; % Oil reserves (billion barrels).
%Plant.Lo = 1.53E3; % Oil reserves (billion barrels).
Plant.prodcapo = 100; % Full annual oil production capacity (billion barrels). Used for the bound constraints in fmincon.
Plant.trigcapo = 0.75; % Ratio of total oil used over oil reserves trigger to reduce oil production capacity as oil is used up.
Plant.Loa = 10; % Amount of resource that can be added annually. Used for the bound constraints in fmincon.
Plant.Loca = 10; % Amount of resource that can be added annually. Used for the bound constraints in fmincon.
%Plant.Loa = 1.53E3; % Annual oil use limit. Used for the bound constraints in fmincon.
Plant.lo = 1; % lead time for new oil field development(years)
Plant.lambdao = 20; % oil field lifetime
lambdao = Plant.lambdao;
Plant.so = 3.56E14; % CO2 emissions per billion barrels of oil used (grams CO2/billion barrels)
%Plant.som = 3.48E11; % methane emissions per billion barrels of oil used (grams CH4/billion barrels)
Plant.som = 0; % methane emissions per billion barrels of oil used (grams CH4/billion barrels)
%Plant.tzo = 31.025; % global oil production capacity at t=0 (billion barrels).
Plant.tzo = 23.27; % portion of global oil production capacity at t=0 allocated to transportation (billion barrels).
Plant.tzoc = 5.0; % portion of global oil production capacity at t=0 allocated to combustion (billion barrels).
mino = 15.0; % Allowable annual decrease in oil production for transportation (billion barrels).
minoc = 6.5; % Allowable annual decrease in oil production for combustion (billion barrels).
AnnCapo = Ko * ( (dis_rate/(((1+dis_rate)^lambdao)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
% Natural Gas Parameters (trillion cubic feet)
Plant.Cg = 5.62E-3*ones(time,1); % variable cost of natural gas ($trillion/ trillion cubic feet)
%Plant.Cg = xlsread('gasvarcost.xls'); % variable cost of natural gas ($trillion/ trillion cubic feet)
Plant.Cge = 5.62E-3*ones(time,1); % variable cost of natural gas for electricity ($trillion/ trillion cubic feet)
%Plant.Cge = xlsread('gaselecvarcost.xls'); % variable cost of natural gas for electricity ($trillion/ trillion cubic feet)
Cg = Plant.Cg;
Cge = Plant.Cge;
%Plant.Kg = 0*ones(time,1); % capital cost of natural gas ($trillion/ trillion cubic feet)
Plant.Kg = 6E-3*ones(time,1); % capital cost of natural gas ($trillion/ trillion cubic feet)
%Plant.Kg = xlsread('gascapcost.xls'); % capital cost of natural gas ($trillion/ trillion cubic feet)
Kg = Plant.Kg;
%Plant.Kge = 1.45E-2*ones(time,1); % capital cost of natural gas electric power plant ($trillion/ trillion cubic feet)
Plant.Kge = 2.05E-2*ones(time,1); % capital cost of natural gas electric power plant ($trillion/ trillion cubic feet)
%Plant.Kge = xlsread('gaseleccapcost.xls'); % capital cost of natural gas electric power plant ($trillion/ trillion cubic feet)
Kge = Plant.Kge;
pg = 0*ones(time,1); % penalty cost of natural gas ($trillion/ trillion cubic feet)
Plant.kg = 1.06; % conversion of natural gas to EJ (EJ/ trillion cubic feet)
Plant.etage = 0.45; % efficiency of burning natural gas to produce electricity.
Plant.etag = 0.77; % efficiency of burning natural gas for specific industrial combustion/furnace needes.
%Plant.Lg = 338; % natural gas reserves (trillion cubic feet).
Plant.Lg = 2.284E4; % natural gas reserves (trillion cubic feet).
%Plant.Lg = 6.84E3; % natural gas reserves (trillion cubic feet).
Plant.prodcapg = 350; % Full annual gas production capacity (trillion cubic feet). Used for the bound constraints in fmincon.
Plant.trigcapg = 0.75; % Ratio of total gas used over gas reserves trigger to reduce gas production capacity as gas is used up.
%Plant.Lga = 6.84E3; % Annual natural gas use limit. Used for the bound constraints in fmincon.
Plant.Lga = 10; % Amount of resource that can be added annually. Used for the bound constraints in fmincon.
Plant.Lgea = 10; % Amount of resource that can be added annually. Used for the bound constraints in fmincon.
Plant.lg = 1; % lead time for new gas field development(years)
Plant.lge = 1; % lead time for new combined cycle electric plants (years).
%Plant.lambdag = 40; % gas field lifetime
Plant.lambdag = 50; % gas field lifetime
lambdag = Plant.lambdag;
Plant.sg = 5.35E13; % CO2 emissions per trillion cubic feet of gas used (grams CO2/trillion cubic feet)
%Plant.sgm = 5.31E11; % methane emissions per trillion cubic feet of gas used (grams CH4/trillion cubic feet)
Plant.sgm = 0; % methane emissions per trillion cubic feet of gas used (grams CH4/trillion cubic feet)
%Plant.tzg = 104.6; % natural gas production capacity at t=0 (trillion cubic feet).
Plant.tzg = 40; % portion of global gas production capacity at t=0 allocated to combustion (trillion cubic feet).
Plant.tzge = 40; % portion of global gas production capacity at t=0 allocated to electricity production (trillion cubic feet).
ming = 1; % Allowable annual increase in gas usage for combustion (trillion cubic feet).
minge = 0.5; % Allowable annual increase in gas usage for electricity (trillion cubic feet).
AnnCapg = Kg * ( (dis_rate/(((1+dis_rate)^lambdag)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
AnnCapge = Kge * ( (dis_rate/(((1+dis_rate)^lambdag)- 1)) + dis_rate ); % Conversion of power plant capital cost to an annual value.
% Coal Parameters (million metric tons)
Plant.Cc = 4.39E-5*ones(time,1); % variable cost of coal ($trillion/ million metric tons)
%Plant.Cc = xlsread('coalvarcost.xls'); % variable cost of coal ($trillion/ million metric tons)
Plant.Cce = 4.39E-5*ones(time,1); % variable cost of coal for electricity ($trillion/ million metric tons)
%Plant.Cce = xlsread('coalelecvarcost.xls'); % variable cost of coal for electricity ($trillion/ million metric tons)
Cc = Plant.Cc;
Cce = Plant.Cce;
Plant.Cct = 6.07E-5*ones(time,1); % variable cost of coal converted to liquid transportation fuel ($trillion/ million metric tons)
%Plant.Cct = xlsread('coalliqvarcost.xls'); % variable cost of coal converted to liquid transportation fuel ($trillion/ million metric tons)
Cct = Plant.Cct;
%Plant.Kc = 0*ones(time,1); % capital cost of coal ($trillion/ million metric tons)
Plant.Kc = 2.37E-2*ones(time,1); % capital cost of coal ($trillion/ million metric tons)
%Plant.Kc = xlsread('coalcapcost.xls'); % capital cost of coal ($trillion/ million metric tons)
Kc = Plant.Kc;
Plant.Kce = 2.46E-2*ones(time,1); % capital cost of coal fired power plant($trillion/ million metric tons)
%Plant.Kce = 9.225E-4*ones(time,1); % capital cost of coal fired power plant($trillion/ million metric tons)
%Plant.Kce = xlsread('coaleleccapcost.xls'); % capital cost of coal fired power plant($trillion/ million metric tons)
Kce = Plant.Kce;
Plant.Kct = 1.99E-1*ones(time,1); % capital cost of coal to liquids conversion ($trillion/ million metric tons)
%Plant.Kct = xlsread('coalliqcapcost.xls'); % capital cost of coal to liquids conversion ($trillion/ million metric tons)
Kct = Plant.Kct;
pc = 0*ones(time,1); % penalty cost of coal ($trillion/ million metric tons)
Plant.kc = 2.93E-2; % conversion of coal to EJ (EJ/ million metric tons)
Plant.etace = 0.36; % efficiency of burning coal to produce electrcity.
Plant.etac = 0.79; % efficiency of burning coal specific industrial combustion/furnace needes.
Plant.etact = 0.16; % efficiency of converting coal to liquid fuel and using it for transportation fuel.
Plant.Lc = 8.62E5; % coal reserves (million metric tons).
Plant.prodcapc = 8000; % Full annual coal production capacity (million metric tons). Used for the bound constraints in fmincon.
Plant.trigcapc = 0.75; % Ratio of total coal used over coal reserves trigger to reduce coal production capacity as coal is used up.
Plant.Lca = 300; % Amount of resource that can be added annually. Used for the bound constraints in fmincon.
Plant.Lcea = 600; % Amount of resource that can be added annually. Used for the bound constraints in fmincon.
Plant.Lcta = 2000; % Amount of resource that can be added annually. Used for the bound constraints in fmincon.
Plant.lc = 1; % lead time for coal development(years)
Plant.lct = 5; % lead time for coal to liquids plant development (years)
%Plant.lambdac = 30; % coal mine, power plant, and coal to liquids plant lifetime
Plant.lambdac = 50; % coal mine, power plant, and coal to liquids plant lifetime
lambdac = Plant.lambdac;
Plant.sc = 2.05E12; % CO2 emissions per million metric tons of coal used (grams CO2/million metric tons)
%Plant.scm = 4.04E9; % methane emissions per million metric tons of coal used (grams CH4/million metric tons)
Plant.scm = 0; % methane emissions per million metric tons of coal used (grams CH4/million metric tons)
Plant.sct = 1.34E12; % Full "well to wheel" CO2-eq emissions associated with coal to liquids for transport fuel (grams CO2-eq/million metric tons)
%Plant.tzc = 4750; % coal production capacity at t=0 (million metric tons).
Plant.tzc = 1615; % coal production capacity at t=0 allocated to combustion(million metric tons).
Plant.tzce = 3088; % coal production capacity at t=0 allocated to electricity production (million metric tons).
Plant.tzct = 0; % coal production capacity at t=0 allocated to coal-to-liquids (million metric tons).
minc = 5; % allowable annual decrease in coal usage for combustion (million metric tons)
mince = 5; % allowable annual decrease in coal usage for electricity (million metric tons)
minct = 4; % allowable annual decrease in coal usage for coal-to-liquids (million metric tons)
Plant.Zc = 4.53E-4*ones(time,1); % decommissioning costs of a coal plant($trillion/million metric tons)
Zc = Plant.Zc;
AnnDecomc = Zc * (dis_rate/(((1+dis_rate)^lambdac)- 1)) ; % Conversion of decommissioning cost to an annual value.
AnnCapc = Kc * ( (dis_rate/(((1+dis_rate)^lambdac)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
AnnCapce = (Kce * ( (dis_rate/(((1+dis_rate)^lambdac)- 1)) + dis_rate) ) + AnnDecomc; % Conversion of power plant capital cost to an annual value.
AnnCapct = Kct * ( (dis_rate/(((1+dis_rate)^lambdac)- 1)) + dis_rate ); % Conversion of coal to liquids plant capital cost to an annual value.
% Nuclear Parameters (Terrawatts)
Plant.Cn = 0.263*ones(time,1); % variable cost of nuclear ($trillion/ TW-year)
%Plant.Cn = xlsread('nukevarcost.xls'); % variable cost of nuclear ($trillion/ TW-year)
Cn = Plant.Cn;
Plant.Kn = 2.3*ones(time,1); % capital cost of nuclear ($trillion/ Terrawatt)
%Plant.Kn = xlsread('nukecapcost.xls'); % capital cost of nuclear ($trillion/ Terrawatt)
Kn = Plant.Kn;
pn = 0*ones(time,1); % penalty cost of nuclear ($trillion/ Terrawatt)
Plant.kn = 31.5; % conversion of nuclear to EJ (EJ/ Terrawatt)
Plant.etan = 0.91; % Nuclear plant capacity factor.
%Plant.etan = 0.60; % Nuclear plant capacity factor.
Plant.Ln = 2.89E3; % Nuclear potential (TW-years of energy). Value based on amount of nuclear fuel available without reprocessing.
%Plant.Lna = 2.89E5; % Annual nuclear use limit. Used for the bound constraints in fmincon.
Plant.Lna = 0.09; % Amount of nuclear that can be added each year (TW). Used to set bound constraints in fmincon.
Plant.ln = 6; % lead time for nuclear plant development(years)
%Plant.lambdan = 40; % nuclear plant lifetime
Plant.lambdan = 50; % nuclear plant lifetime
lambdan = Plant.lambdan;
Plant.sn = 6.5E10; % CO2 emissions per TW-year of nuclear used (grams CO2/TW-year). Only needed when new capacity is added.
Plant.tzn = 3.49E-1; % installed nuclear capacity at t=0 (TW).
minn = 1E-5; % allowable annual decrease in nuclear usage (TW).
Plant.Zn = 3.5E-1*ones(time,1); % decommissioning costs ($trillion/Terrawat)
Zn = Plant.Zn;
AnnDecomn = Zn * (dis_rate/(((1+dis_rate)^lambdan)- 1)) ; % Conversion of decommissioning cost to an annual value.
AnnCapn = (Kn * ( (dis_rate/(((1+dis_rate)^lambdan)- 1)) + dis_rate )) + AnnDecomn; % Conversion of capital + decomm cost to an annual value.
% Hydro Parameters (Terrawatts)
Plant.Ch = 0*ones(time,1); % variable cost of hydro ($trillion/ TW-year)
%Plant.Ch = xlsread('hydrovarcost.xls'); % variable cost of hydro ($trillion/ TW-year)
Ch = Plant.Ch;
Plant.Kh = 2.95*ones(time,1); % capital cost of hydro ($trillion/ TW-year)
%Plant.Kh = xlsread('hydrocapcost.xls'); % capital cost of hydro ($trillion/ TW-year)
Kh = Plant.Kh;
ph = 0*ones(time,1); % penalty cost of hydro ($trillion/ TW-year)
Plant.kh = 31.5; % conversion of terrawatt-y of hydro to EJ (EJ/TW-y)
Plant.etah = 0.40; % hydro plant capacity factor
%Plant.etah = 0.30; % hydro plant capacity factor
%Plant.Lh = 0.157; % global hydro power potential (TW).
Plant.Lh = 3.72; % global hydro power potential (TW).
%Plant.Lha = 3.72; % Annual hydro use limit. Used for the bound constraints in fmincon.
Plant.Lha = 0.02; % Amount of resource that can be added annually. Use to compute annual hydro use limit for the bound constraints in fmincon.
Plant.lh = 4; % lead time for hydro plant development(years)
Plant.lambdah = 75; % hydro plant lifetime
lambdah = Plant.lambdah;
Plant.sh = 1E10; % CO2 emissions per TW of hydro used (grams CO2/TW). Only needed when new capacity is added.
Plant.tzh = 0.7; % installed hydro capacity at t=0 (TW).
minh = 1E-4; % allowable annual decrease in hydro usage (TW).
Plant.Zh = 8.85E-1*ones(time,1); % decommissioning cost of hydro ($trillion/TW)
Zh = Plant.Zh;
AnnDecomh = Zh * (dis_rate/(((1+dis_rate)^lambdah)- 1)) ; % Conversion of decommissioning cost to an annual value.
AnnCaph = (Kh * ( (dis_rate/(((1+dis_rate)^lambdah)- 1)) + dis_rate )) + AnnDecomh; % Conversion of capital + decomm cost to an annual value.
% Wind Parameters (Terrawatts)
Plant.Cw = 0*ones(time,1); % variable cost of wind ($trillion/ TW-year)
%Plant.Cw = xlsread('windvarcost.xls'); % variable cost of wind ($trillion/ TW-year)
Cw = Plant.Cw;
Plant.Kw = 2.25*ones(time,1); % capital cost of wind ($trillion/ Terrawatt)
%Plant.Kw = 2.0*ones(time,1); % capital cost of wind ($trillion/ Terrawatt) (assume a subsidy)
%Plant.Kw = xlsread('windcapcost.xls'); % capital cost of wind ($trillion/ Terrawatt)
Kw = Plant.Kw;
pw = 0*ones(time,1); % penalty cost of wind ($trillion/ Terrawatt)
Plant.kw = 31.5; % conversion of wind to EJ (EJ/ Terrawatt)
Plant.etaw = 0.30; % wind plant capacity factor.
%Plant.etaw = 0.25; % wind plant capacity factor.
%Plant.Lw = 11; % wind potential (TW).
Plant.Lw = 196; % wind potential (TW).
%Plant.Lwa = 196; % Annual wind use limit. Used for the bound constraints in fmincon.
Plant.Lwa = 0.1; % Amount of resource that can be added annually. Use to compute annual wind use limit for the bound constraints in fmincon.
%Plant.Lwa = 0.55; % Amount of resource that can be added annually. Use to compute annual wind use limit for the bound constraints in fmincon.
Plant.lw = 1; % lead time for wind plant development(years)
Plant.lambdaw = 25; % wind plant lifetime
lambdaw = Plant.lambdaw;
Plant.sw = 1.2E10; % CO2 emissions per TW of wind used (grams CO2/TW). Only needed when new capacity is added.
Plant.tzw = 1.74E-2; % global installed wind capacity at t=0 (TW).
%Plant.tzw = 6.68E-2; % installed wind capacity at t=0 (TW).
minw = 1.0E-4; % allowable annual decrease in wind usage (TW).
AnnCapw = Kw * ( (dis_rate/(((1+dis_rate)^lambdaw)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
% Geothermal Parameters (Terrawatts)
Plant.Cgt = 0*ones(time,1); % variable cost of geothermal ($trillion/ TW-year)
%Plant.Cgt = xlsread('geothermvarcost.xls'); % variable cost of geothermal ($trillion/ TW-year)
Cgt = Plant.Cgt;
Plant.Kgt = 4.46*ones(time,1); % capital cost of geothermal ($trillion/ Terrawatt)
%Plant.Kgt = xlsread('geothermcapcost.xls'); % capital cost of geothermal ($trillion/ Terrawatt)
Kgt = Plant.Kgt;
pgt = 0*ones(time,1); % penalty cost of geothermal ($trillion/ Terrawatt)
Plant.kgt = 31.5; % conversion of geothermal to EJ (EJ/ Terrawatt)
Plant.etagt = 0.81; % geothermal plant capacity factor.
%Plant.Lgt = 2E-2; % geothermal potential (TW).
Plant.Lgt = 0.209; % geothermal potential (TW).
%Plant.Lgta = 0.209; % Annual geothermal use limit. Used for the bound constraints in fmincon.
Plant.Lgta = 0.004; % Amount of resource that can be added annually. Use to compute annual geothermal use limit for the bound constraints in fmincon.
Plant.lgt = 4; % lead time for geothermal plant development(years)
Plant.lambdagt = 25; % geothermal plant lifetime
lambdagt = Plant.lambdagt;
Plant.sgt = 9.1E10; % CO2 emissions per TW of geothermal used (grams CO2/TW). Only needed when new capacity is added.
Plant.tzgt = 7.97E-3; % installed geothermal capacity at t=0 (TW).
mingt = 1E-5; % allowable annual decrease in geothermal usage (TW).
AnnCapgt = Kgt * ( (dis_rate/(((1+dis_rate)^lambdagt)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
% Solar-PV Parameters (Terrawatts)
Plant.Csp = 0*ones(time,1); % variable cost of solar-PV ($trillion/ TW-year)
%Plant.Csp = xlsread('solarpvvarcost.xls'); % variable cost of solar-PV ($trillion/ TW-year)
Csp = Plant.Csp;
Plant.Ksp = 4.21*ones(time,1); % capital cost of solar-PV ($trillion/ Terrawatt)
%Plant.Ksp = 2.5*ones(time,1); % capital cost of solar-PV ($trillion/ Terrawatt) (assume a subsidy)
%Plant.Ksp = xlsread('solarpvcapcost.xls'); % capital cost of solar-PV ($trillion/ Terrawatt)
Ksp = Plant.Ksp;
psp = 0*ones(time,1); % penalty cost of solar-PV ($trillion/ Terrawatt)
Plant.ksp = 31.5; % conversion of solar-PV to EJ (EJ/ Terrawatt)
Plant.etasp = 0.20; % solar-PV capacity factor.
Plant.Lsp = 2.0; % solar-PV potential (TW).
%Plant.Lspa = 2; % Annual solar-PV use limit. Used for the bound constraints in fmincon.
Plant.Lspa = 5.0E-3; % Amount of resource that can be added annually. Use to compute annual solar-PV use limit for the bound constraints in fmincon.
Plant.lsp = 1; % lead time for solar-PV plant development(years)
Plant.lambdasp = 25; % solar-PV lifetime
lambdasp = Plant.lambdasp;
Plant.ssp = 6E10; % CO2 emissions per TW of solar-PV used (grams CO2/TW). Only needed when new capacity is added.
Plant.tzsp = 1.25E-3; % global installed solar-PV capacity at t=0 (TW).
%Plant.tzsp = 1.55E-2; % installed solar-PV capacity at t=0 (TW).
minsp = 1.0E-5; % allowable annual decrease in solar-PV usage (TW).
AnnCapsp = Ksp * ( (dis_rate/(((1+dis_rate)^lambdasp)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
% Solar-thermal Parameters (Terrawatts)
Plant.Cst = 0*ones(time,1); % variable cost of solar-thermal ($trillion/ TW-year)
% Plant.Cst = xlsread('solarthermvarcost.xls'); % variable cost of solar-thermal ($trillion/ TW-year)
Cst = Plant.Cst;
Plant.Kst = 5.11*ones(time,1); % capital cost of solar-thermal ($trillion/ Terrawatt)
%Plant.Kst = xlsread('solarthermcapcost.xls'); % capital cost of solar-thermal ($trillion/ Terrawatt)
Kst = Plant.Kst;
pst = 0*ones(time,1); % penalty cost of solar-thermal ($trillion/ Terrawatt)
Plant.kst = 31.5; % conversion of solar-thermal to EJ (EJ/ Terrawatt)
Plant.etast = 0.56; % solar-thermal capacity factor.
%Plant.etast = 0.3; % solar-thermal capacity factor.
Plant.Lst = 1.39; % solar-thermal potential (TW).
%Plant.Lsta = 1.39; % Annual solar-thermal use limit. Used for the bound constraints in fmincon.
Plant.Lsta = 2E-3; % Amount of resource that can be added annually. Use to compute annual solar-thermal use limit for the bound constraints in fmincon.
Plant.lst = 3; % lead time for solar-thermal plant development(years)
Plant.lambdast = 40; % solar-thermal lifetime
lambdast = Plant.lambdast;
Plant.sst = 9E10; % CO2 emissions per TW of solar-thermal used (CO2/TW). Only needed when new capacity is added.
Plant.tzst = 3.54E-4; % global installed solar-thermal capacity at t=0 (TW).
%Plant.tzst = 1.44E-3; % installed solar-thermal capacity at t=0 (TW).
minst = 1.0E-6; % allowable annual decrease in solar-thermal usage (TW).
AnnCapst = Kst * ( (dis_rate/(((1+dis_rate)^lambdast)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
% Biofuel Parameters (Trillion Liters)
Plant.Cbf = 6.57E-1*ones(time,1); % variable cost of biofuel ($trillion/ trillion liters)
%Plant.Cbf = xlsread('biofuelvarcost.xls'); % variable cost of biofuel ($trillion/ trillion liters)
Cbf = Plant.Cbf;
Plant.Kbf = 0.536*ones(time,1); % capital cost of biofuel ($trillion/ trillion liters)
%Plant.Kbf = xlsread('biofuelcapcost.xls'); % capital cost of biofuel ($trillion/ trillion liters)
Kbf = Plant.Kbf;
pbf = 0*ones(time,1); % penalty cost of biofuel ($trillion/trillion liters)
Plant.kbf = 34; % conversion of biofuel to EJ (EJ/trillion liters)
Plant.etabf = 0.23; % efficiency of burning biofuel for transportation.
Plant.etabfc = 0.67; % efficiency of burning biofuel for specific industrial combustion/furnace needs.
Plant.Lbf = 3E3; % acreage available for biofuel production (million ha).
Plant.Lbfa = 3.2; % Amount of resource that can be added annually. Use to compute annual biofuel limit for the bound constraints in fmincon.
Plant.lbf = 1; % lead time for biofuel(years)
Plant.lambdabf = 30; % biofuel plant lifetime
lambdabf = Plant.lambdabf;
Plant.aebf = 2.81E2; % area-energy conversion factor (million ha/trillion liters)
Plant.sbf = 1.40E15; % CO2 emissions per trillion liters of biofuel (grams CO2/trillion liters)
Plant.sbfq = 7.0E14; % CO2 sequestered per trillion liters of biofuel (grams CO2/trillion liters). Only used in reserves calculation
Plant.tzbf = 1.81E-2; % biofuel production capacity at t=0 (trillion liters).
Plant.tzbfc = 2E-4; % global biofuel production capacity at t=0 allocated to combustion (trillion liters).
minbf = 1.0E-3; % allowable annual decrease in biofuel usage for transportation (trillion liters).
minbfc = 1E-4; % allowable annual decrease in biofuel usage for combustion (trillion liters)
AnnCapbf = Kbf * ( (dis_rate/(((1+dis_rate)^lambdabf)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
% Biomass Parameters (Terrawatts)
Plant.Cbm = 5.26E-3*ones(time,1); % variable cost of biomass ($trillion/TW)
%Plant.Cbm = xlsread('biomassvarcost.xls'); % variable cost of biomass ($trillion/TW)
Cbm = Plant.Cbm;
Plant.Kbm = 4.22*ones(time,1); % capital cost of biomass ($trillion/TW)
%Plant.Kbm = xlsread('biomasscapcost.xls'); % capital cost of biomass ($trillion/TW)
Kbm = Plant.Kbm;
pbm = 0*ones(time,1); % penalty cost of biomass ($trillion/TW)
Plant.kbm = 31.5; % conversion of biomass to EJ (EJ/TW)
Plant.etabm = 0.83; % Biomass plant capacity factor
%Plant.etabm = 0.7; % Biomass plant capacity factor
Plant.Lbm = 3E3; % acreage available for biomass production (million ha).
Plant.Lbma = 1E-3; % Amount of resource that can be added annually. Use to compute annual biomass limit for the bound constraints in fmincon.
Plant.lbm = 1; % lead time for biomass plant (years)
Plant.lambdabm = 20; % biomass plant lifetime
lambdabm = Plant.lambdabm;
Plant.aebm = 1.73E2; % area-energy conversion factor (million ha/TW)
Plant.sbm = 1.58E14; % CO2 emissions per TW of biomass (grams CO2/TW)
Plant.sbmq = 1.58E14; % CO2 sequestered per TW of biomass (grams CO2/TW). Only used in reserves calculation
Plant.tzbm = 3.58E-3; % global biofuel production capacity at t=0 (TW).
minbm = 1.0E-3; % allowable annual decrease in biomass usage (TW).
AnnCapbm = Kbm * ( (dis_rate/(((1+dis_rate)^lambdabm)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
% Energy Efficiency Parameters (EJ)
Plant.GDP = 15.7; % Starting year GDP. Use for energy efficiency cost computation. (trillion dollars)
Plant.Ce = 0*ones(time,1); % Placeholder for variable costs.
Ce = Plant.Ce;
Plant.Ke = 1.57E-5*ones(time,1); % Cost of energy efficiency as a fraction of GDP. Must multiply this value value by GDP. ($trillion/EJ)
%Plant.Ke = xlsread('effcapcost.xls'); % Cost of energy efficiency as a fraction of GDP. Must multiply this value value by GDP. ($trillion/EJ)
Ke = Plant.Ke;
pe = 0*ones(time,1); % Placeholder for penalty costs.
%Plant.Le = 0.1*ones(time,1); % Total available energy efficiency as a fraction of total global energy demand.
Plant.Le = xlsread('effamount.xls'); % Total available energy efficiency as a fraction of total global energy demand.
Plant.le = 1; % Lead time for energy efficiency.
Plant.lambdae = 15; % Average lifetime of efficiency measures before they need to be retrofitted.
lambdae = Plant.lambdae;
Plant.tze = 1E-6; % global energy efficiency at t=0.
mine = 10; % allowable annual energy efficiency usage (EJ)
AnnCape = Ke * ( (dis_rate/(((1+dis_rate)^lambdae)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
% Sulfur Injection Parameters (Million metric tons SO2)
Plant.Cs = 4.55E-4*ones(time,1); % Variable cost of sulfur injection ($trillion/million metric tons).
%Plant.Cs = xlsread('sulfurvarcost.xls'); % Variable cost of sulfur injection ($trillion/million metric tons).
Cs = Plant.Cs;
Plant.Ks = 1.05E-3*ones(time,1); % Capital cost of sulfur injection ($trillion/million metric tons).
%Plant.Ks = xlsread('sulfurcapcost.xls'); % Capital cost of sulfur injection ($trillion/million metric tons).
Ks = Plant.Ks;
ps = 0*ones(time,1); % Penalty cost of sulfur injection
Plant.Ls = 20; % Limit in million metric tons based on number of planes needed to replicate Mt. Pinatubo sulfur injection
Plant.Lsa = 2; % Allowable annual increase in sulfur injection deployment.
Plant.ls = 1; % lead time for sulfur injection (years)
Plant.lambdas = 3; % time sulfur injection remains effective (years). Use in climate equation
Plant.lambdas2 = 50; % lifetime of air tankers needed to deploy the sulfur. Use this in the cost equation.
lambdas2 = Plant.lambdas2;
Plant.ss = 0.75; % Radiative forcing reduction due to sulfur injection (W/sq.meter/million metric tons)
Plant.tzs = 1E-6; % global sulfur injection capacity at t=0.
mins = 1E6; % allowable annual decrease in sulfur usage.
AnnCaps = Ks * ( (dis_rate/(((1+dis_rate)^lambdas2)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
% Iron Seeding Parameters (million Metric Tons)
Plant.Ci = 6.75E-3*ones(time,1); % Variable cost of iron seeding ($trillion/million metric tons).
%Plant.Ci = xlsread('ironvarcost.xls'); % Variable cost of iron seeding ($trillion/million metric tons).
Ci = Plant.Ci;
Plant.Ki = 5.8E-4*ones(time,1); % Capital cost of new ships for iron seeding ($trillion/million metric tons)
%Plant.Ki = xlsread('ironcapcost.xls'); % Capital cost of new ships for iron seeding ($trillion/million metric tons)
Ki = Plant.Ki;
pi = 0*ones(time,1); % Penalty cost of iron seeding.
Plant.Li = 0.6; % Limit on the amount of feasible iron seeding (million metric tons).
Plant.Lia = 0.05; % Allowable annual increase in iron seeding deployment.
Plant.li = 1; % lead time for iron seeding (years)
Plant.lambdai = 1; % time iron seeding remains effective (years). Use in climate equation.
Plant.lambdai2 = 30; % life time of tankers (years). Use in cost equations
lambdai2 = Plant.lambdai2;
Plant.si = 1.02E16; % CO2 sequestered per million metric tons of iron seeding (grams CO2/milllion metric tons)
Plant.tzi = 1E-6; % global iron seeding capacity at t=0.
mini = 0.7; % allowable decrease in iron seeding usage.
AnnCapi = Ki * ( (dis_rate/(((1+dis_rate)^lambdai2)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
% Tree Planting Parameters (million ha)
Plant.Ct = 0*ones(time,1); % Placeholder for variable costs.
Ct = Plant.Ct;
Plant.Kt = 1.50E-3*ones(time,1); % Capital cost of tree planting ($trillion/million ha)
%Plant.Kt = xlsread('treescapcost.xls'); % Capital cost of tree planting ($trillion/million ha)
Kt = Plant.Kt;
pt = 0*ones(time,1); % Placeholder for penalty costs
Plant.Lt = 3E3; % Acreage available for tree planting (million ha)
Plant.Lta = 1; % Limit on acreage of trees that can be planted annually. Use to compute annual tree limit for the bound constraints in fmincon.
Plant.lt = 1; % Lead time for tree planting to take effect
Plant.steep = 1; % Steepness parameter for the logistic function .
steep = Plant.steep;
Plant.midpoint = 25; % Midpoint of the logistic function.
midpoint = Plant.midpoint;
dlogistic = exp((midpoint-1)*steep)/((1 + exp((midpoint-1)*steep))^2); % derivative of the logistic function at t = 1.
Plant.dlogistic = dlogistic;
Plant.lambdat = 70; % time that tree planting stays effective (years)
lambdat = Plant.lambdat;
Plant.st = 2.6E13; % CO2 sequestered per million ha (grams CO2/milllion ha)
Plant.tzt = 50; % global forest area at t=0 (million ha).
mint = 3E3; % allowable annual reduction in tree planting.
AnnCapt = Kt * ( (dis_rate/(((1+dis_rate)^lambdat)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
% Sea Spray Cloud Enhancement Parameters (cu. meter/sec)
Plant.Css = 0*ones(time,1); % Variable cost of sea spray. Needed as a place-holder to form the matrices below
Css = Plant.Css;
Plant.Kss = 1.04E-4*ones(time,1); % Capital cost of sea spray geoengineering ($trillion/cu. meter per sec)
%Plant.Kss = xlsread('seaspraycapcost.xls'); % Capital cost of sea spray geoengineering ($trillion/cu. meter per sec)
Kss = Plant.Kss;
pss = 0*ones(time,1); % Penalty cost of sea spray
Plant.Lss = 1; % Limit based on declining effectiveness of sea spray with increased rate (in cu. meter/sec).
% Max sea spray limit also set due to its high equivalent emissions value that would otherwise cause bad
% results.
Lss = Plant.Lss;
Plant.Lssa = 10; % Allowable annual increase in sea spray deployment
Plant.lss = 1; % lead time for sea spray based on time to construct ships (years)
Plant.lambdass = 1; % time sea spray remains effective (years). Use in climate equation
Plant.lambdass2 = 10; % lifetime of ships needed to deploy the sea spray. Use this in the cost equation.
lambdass2 = Plant.lambdass2;
Plant.ssf = 8.22E-2; % Radiative forcing reduction due to sea spray (W/sq.meter/cu. meter per sec)
Plant.Ess = 9.65E16; % Equivalent emissions reduction (grams CO2/cu meter per sec).
Plant.tzss = 1E-6; % global sea spray capacity at t=0.
minss = 100; % allowable annual decrease in sea spray usage.
AnnCapss = Kss * ( (dis_rate/(((1+dis_rate)^lambdass2)- 1)) + dis_rate ); % Conversion of capital cost to an annual value.
% create vectors of all the variable, penalty, capital, and annual capital costs, the combined costs, & the lambda values.
VarCosts = [Co Cg Cc Cn Ch Cw Cgt Csp Cst Cbf Cbm Ce Cs Ci Ct Css Co Cge Cce Cct Cbf];
Plant.VarCosts = VarCosts;
penalty = [po pg pc pn ph pw pgt psp pst pbf pbm pe ps pi pt pss po pg pc pc pbf];
Plant.penalty = penalty;
AnnCosts = [AnnCapo AnnCapg AnnCapc AnnCapn AnnCaph AnnCapw AnnCapgt AnnCapsp AnnCapst AnnCapbf AnnCapbm AnnCape AnnCaps AnnCapi AnnCapt AnnCapss...
AnnCapo AnnCapge AnnCapce AnnCapct AnnCapbf];
Plant.AnnCosts = AnnCosts;
CapCosts = [Ko Kg Kc Kn Kh Kw Kgt Ksp Kst Kbf Kbm Ke Ks Ki Kt Kss...
Ko Kge Kce Kct Kbf];
Plant.CapCosts = CapCosts;
Plant.T = VarCosts + AnnCosts + penalty;
% Total of variable plus annualized capital and decommissioning costs.
Plant.lambda = [lambdao lambdag lambdac lambdan lambdah lambdaw lambdagt lambdasp lambdast lambdabf lambdabm lambdae lambdas2 lambdai2 lambdat lambdass2...
lambdao lambdag lambdac lambdac lambdabf];
% Matrix of adjustments for any unrecovered annualized capital costs
% due to resource usage peaking, then dropping off.
adjo = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/billion barrels).
adjg = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/trillion cubic feet).
adjc = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/million metric tons).
adjn = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/TW).
adjh = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/TW).
adjw = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/TW).
adjgt = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/TW).
adjsp = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/TW).
adjst = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/TW).
adjbf = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/trillion liters).
adjbm = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/TW).
adje = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/EJ).
adjs = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/million tons).
adji = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/million metric tons).
adjt = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/million ha).
adjss = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/cu.meter/sec).
adjoc = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/billion barrels).
adjge = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/trillion cubic feet).
adjce = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/million metric tons).
adjct = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/million metric tons).
adjbfc = 0*ones(time,1); % adjustment cost for unrecovered cost due to resource peak ($trillion/trillion liters).
Plant.adj = [adjo adjg adjc adjn adjh adjw adjgt adjsp adjst adjbf adjbm adje adjs adji adjt adjss...
adjoc adjge adjce adjct adjbfc];
% Specify the capacities (not the actual values) of the energy and geoengineering resources at time t=0 in order to
% determine what new incremental capacities are added in year 1 and to
% insure that lead times for new resource additions are respected.
tzo = Plant.tzo;
tzoc = Plant.tzoc;
tzg = Plant.tzg;
tzge = Plant.tzge;
tzc = Plant.tzc;
tzce = Plant.tzce;
tzct = Plant.tzct;
tzn = Plant.tzn;
tzh = Plant.tzh;
tzw = Plant.tzw;
tzgt = Plant.tzgt;
tzsp = Plant.tzsp;
tzst = Plant.tzst;
tzbf = Plant.tzbf;
tzbfc = Plant.tzbfc;
tzbm = Plant.tzbm;
tze = Plant.tze;
tzs = Plant.tzs;
tzi = Plant.tzi;
tzt = Plant.tzt;
tzss = Plant.tzss;
Plant.tzero = [tzo tzg tzc tzn tzh tzw tzgt tzsp tzst tzbf tzbm tze tzs tzi tzt...
(tzo+tzoc) (tzg+tzge) (tzc+tzce+tzct) tzh tzgt (tzbf+tzbfc) tzbm...
tzss tzoc tzge tzce tzct tzbfc];
% Create matrix of initial values of the non-reserve resources used to
% calculate the lower bound LoB at each time step.
Plant.initial = [tzo tzg tzc tzn tzh tzw tzgt tzsp tzst tzbf tzbm tze tzs tzi tzt...
tzss tzoc tzge tzce tzct tzbfc];
% Create matrix of allowable annual resource reductions.
Plant.min = [mino ming minc minn minh minw mingt minsp minst minbf minbm mine mins mini mint minss...
minoc minge mince minct minbfc];
% Specify the initial values of oil(transport), natural gas(combustion), coal(combustion), nuclear, hydro, wind,
% geothermal, solar-PV, solar-thermal, biofuel(transport), biomass, energy efficiency,
% sulfur injection, iron seeding, tree planting, reserve oil, reserve gas,
% reserve coal, reserve hydro, reserve geothermal, reserve biofuel,
% reserve biomass, oil(combustion), gas(electric), coal(electric),
% coal(transport), and biofuel(combustion)(in that order) for use in the optimization algorithm.
R = Plant.R;
Plant.StartingValue = [tzo tzg tzc tzn tzh tzw tzgt tzsp tzst tzbf tzbm tze tzs tzi tzt...
(R*(tzo+tzoc)) (R*(tzg+tzge)) (R*(tzc+tzce+tzct)) (R*tzh) (R*tzgt) (R*(tzbf+tzbfc)) (R*tzbm) tzss tzoc tzge...
tzce tzct tzbfc];
% set up vectors of each resource and geoengineering option so that their values can be
% stored at the end of each time step and be available for later
% calculations. Since we are looking 100 years, each vector will
% contain 100 entries. Also set up a vector for keeping track of the
% the computed energy.
Plant.qo = zeros(time,1);
Plant.qg = zeros(time,1);
Plant.qc = zeros(time,1);
Plant.qn = zeros(time,1);
Plant.qh = zeros(time,1);
Plant.qw = zeros(time,1);
Plant.qgt = zeros(time,1);
Plant.qsp = zeros(time,1);
Plant.qst = zeros(time,1);
Plant.qbf = zeros(time,1);
Plant.qbm = zeros(time,1);
Plant.qe = zeros(time,1);
Plant.qs = zeros(time,1);
Plant.qi = zeros(time,1);
Plant.qt = zeros(time,1);
Plant.qor = zeros(time,1); % oil held as energy reserves
Plant.qgr = zeros(time,1); % natural gas held as energy reserves
Plant.qcr = zeros(time,1); % coal held as energy reserves
Plant.qhr = zeros(time,1); % hydro held as energy reserves
Plant.qgtr = zeros(time,1); % geothermal held as energy reserves
Plant.qbfr = zeros(time,1); % biofuel held as energy reserves
Plant.qbmr = zeros(time,1); % biomass held as energy reserves
Plant.qss = zeros(time,1);
Plant.qoc = zeros(time,1);
Plant.qge = zeros(time,1);
Plant.qce = zeros(time,1);
Plant.qct = zeros(time,1);
Plant.qbfc = zeros(time,1);
% Form full array of the final resource values.
Plant.q = zeros(time,TotalResources);
% Set up matrix to keep track of how much of each resource was
% incrementally added at each time step (both due to increased resource need and
% resources reaching the end of their service lives). Also set up a
% matrix to track of how much resource is retiring (including resources
% that were used for reserves).
Plant.qdot = zeros(time,Resources); % change in resources combines changes due to deployment and reserves.
Plant.retire = zeros(time,Resources); % amount of resources reaching the end of their service lives.
% Set up matrix of final resource values with the resources re-ordered
% so that subcomponents are grouped together. This is just for
% ease of reading the Excel printouts.
Plant.printq = zeros(time,TotalResources);
Plant.Q = ones(time,1); % Total of all energy produced
Plant.Ql = ones(time,1); % Total energy produced from liquid fuels.
Plant.Qe = ones(time,1); % Total energy produced from electricity.
Plant.Qc = ones(time,1); % Total energy produced from combustion fuels.
Plant.Reserves = ones(time,1); % Energy reserves at each time step
% Set up additional vectors to keep track of overall emissions (S), CO2
% concentration (c), radiative forcing (dF, dFco2, dFsulfur, dFseaspray),
% temperature change (dT),
% and total cost. Set the concentration co2
% initially equal to the initial atmospheric CO2 concentration to avoid
% possible divide by zero and square root problems.
Plant.S = zeros(time,1); % Total CO2 emissions in grams
Plant.Sm = zeros(time,1); % Total methane emissions in grams
Plant.Sn = zeros(time,1); % N2O emissions in CO2-eq.
Plant.Stotal = zeros(time,1); % Total CO2-eq emissions, including effects of tree planting and iron seeding.
Plant.Sadjust = zeros(time,1); % Total CO2-eq emissions, including effect of sulfur geoengineering.
Plant.co2 = Plant.co2i * ones(time,1); % CO2 concentrations (ppm)
Plant.ch4 = Plant.ch4i * ones(time,1); % CH4 concentrations (ppb)
Plant.n2o = Plant.n2oi * ones(time,1); % N2O concentrations (ppb)
Plant.co2e = zeros(time,1); % Equivalent CO2 concentration in ppm after including effects of CH4 & any sulfur injection.
Plant.emitdecayCO2 = zeros(time,1); % Matrix of decays of each CO2 emissions pulse
Plant.emitdecayCH4 = zeros(time,1); % Matrix of decays of each CH4 emissions pulse
Plant.emitdecayN2O = zeros(time,1); % Matrix of decays of each N2O emissions pulse
Plant.dco2 = zeros(time,1); % Change in CO2 concentration in ppm
Plant.dch4 = zeros(time,1); % Change in CH4 concentration in ppb
Plant.dn2o = zeros(time,1); % Change in N2O concentration in ppb
Plant.dF = zeros(time,1);
Plant.dFco2 = zeros(time,1);
Plant.dFch4 = zeros(time,1);
Plant.dFn2o = zeros(time,1);
Plant.dFsulfur = zeros(time,1);
Plant.dFseaspray = zeros(time,1);
Plant.dT = zeros(time,1);
Plant.Filename = 'Toy_Problem_Plant'; % designate name for the plant data file
save(Plant.Filename,'Plant');