diff --git a/Snakefile b/Snakefile
index c0cbd90..49954a2 100644
--- a/Snakefile
+++ b/Snakefile
@@ -28,7 +28,8 @@ rule retrieve_load:
 rule retrieve_existing_generators:
     output: 
         generators = "data/existing_generators.csv",
-        aggregated = "data/aggregated_generators.csv"
+        aggregated = "data/aggregated_generators.csv",
+        build_year = "data/build_year_generators.csv"
     script: "scripts/retrieve_generators.py"
 
 rule retrieve_renewable_profiles:
@@ -39,6 +40,11 @@ rule retrieve_renewable_profiles:
         solar = "data/time_series/solar.csv"
     script: "scripts/retrieve_renewables.py"
 
+rule retrieve_co2_emissions:
+    output: 
+        emissions = "data/technology_emissions.csv"
+    script: "scripts/retrieve_emissions.py"
+
 rule build_topology:
     input: 
         supply_regions = "data/spatial_data/supply_regions.shp"
@@ -63,17 +69,27 @@ rule add_electricity:
         costs = "data/technology_costs.csv",
         wind_profile = "data/time_series/wind.csv",
         solar_profile = "data/time_series/solar.csv",
-        base_network = "data/networks/base_network.nc"
+        base_network = "data/networks/base_network.nc",
+        emissions = "data/technology_emissions.csv",
+        build_years = "data/build_year_generators.csv"
     output:
         elec_network = "data/networks/electricity_network.nc",
         final_costs = "data/final_costs.csv"
     script: "scripts/add_electricity.py"
 
-
 rule solve_network:
     input:
         elec_network = "data/networks/electricity_network.nc"
     output:
-        solved_network = "results/networks/illinois_solved.nc",
-        dispatch_figure = "results/figures/illinois_dispatch.png"
-    script: "scripts/solve_network.py"
\ No newline at end of file
+        solved_network = "results/networks/illinois_solved.nc"
+    script: "scripts/solve_network.py"
+
+rule plot_results:
+    input:
+        solved_network = "results/networks/illinois_solved.nc"
+    output: 
+        dispatch_figure = "results/figures/illinois_dispatch.png",
+        capacity_figure = "results/figures/illinois_capacity.png",
+        emissions_figure = "results/figures/illinois_emissions.png",
+        active_units_figure = "results/figures/illinois_active_units.png"
+    script: "scripts/plot_results.py"
\ No newline at end of file
diff --git a/config.yml b/config.yml
index 4c448bb..7b49e92 100644
--- a/config.yml
+++ b/config.yml
@@ -1,12 +1,27 @@
 state_abbr: 'IL'
 
 geo_res: 'rto'  # accepts: 'rto' or 'county'
-time_res: 4  # hours 
+time_res: 1  # hours 
+load_filter: 60e3  # MW, load above this level will be removed as outliers.
+total_demand: 136e6  # Annual MWh in the first year
+load_growth: 0.00  # % annual growth
+load_share:
+  MISO-Z4: 0.33
+  ComEd: 0.67
 
-model_year: 2025
+model_years: [2025, 2030, 2035, 2040, 2045]
+plot_year: 2045
 
-solar_year: 2019  # NSRDB only accepts years 1998-2020
-wind_year: 2011  # WTK only goes from 2009-2013
+co2_limits:
+  2020: 31.25
+  2025: 25  # MT
+  2030: 18.75
+  2035: 12.5
+  2040: 6.25
+  2045: 0.0
+
+solar_years: [2016,2017,2018,2019,2020]  # NSRDB only accepts years 1998-2020
+wind_years: [2009,2010,2011,2012,2013]  # WTK only goes from 2009-2013
 load_year: 2019
 
 rto_subba: ['0004','CE']  # MISO-Z4, ComEd
@@ -14,14 +29,25 @@ region_names: ['MISO-Z4','ComEd']
 
 discount_rate: 0.07
 
+growth_rates:  # MW/year
+  Nuclear: 1e3
+  Biomass: 1e3
+  Natural Gas: 10e3
+  Coal: 10e3
+  Petroleum: 10e3
+  Batteries: 5e3
+  Solar: 5e3
+  Wind: 5e3
+
 lifetime:
   Nuclear: 80
   Biomass: 40
   Natural Gas: 40
-  Coal: 40
+  Coal: 50
   Solar: 20
   Wind: 20
   Batteries: 15
+  Petroleum: 40
 
 atb_params:
     atb_year: 2022  # the ATB publication year
@@ -56,11 +82,10 @@ atb_params:
 lines:  # this is an assumption
   v_nom: 400.
   s_nom: 0.
-  x: 1  # reactance
-  r: 1  # resistance
+  x: 0.001  # reactance
+  r: 0.001  # resistance
   s_nom_extendable: True
 
-
 turbine_params:
   cut_in: 3.0  # m/s
   cut_out: 25.0  # m/s
@@ -69,35 +94,25 @@ turbine_params:
   rated_power: 2.75  # MW
   air_density: 1.225  # kg/m^3
 
-# extendable_techs:  [
-#       'CCAvgCF', 
-#       'CTAvgCF', 
-#       'Biopower', 
-#       'Land-Based Wind',
-#       'Utility PV', 
-#       'IGCCAvgCF', 
-#       'LWR',
-#       'NuclearSMR', 
-#       '10Hr Battery Storage', 
-#       '8Hr Battery Storage',
-#       '6Hr Battery Storage', 
-#       '4Hr Battery Storage', 
-#       '2Hr Battery Storage'
-#       ]
+nuclear_params:  # for existing reactors, 2023 NEI costs in context
+  capacity_factor: 0.972
+  capital_cost: 6.88 # $/MWh
+  fixed_om: 18.68 # $/MWh
+  fuel: 5.37  # $/MWh
 
 extendable_techs:  [
       # 'CCAvgCF', 
       # 'CTAvgCF', 
-      # 'Biopower', 
-      # 'Land-Based Wind',
-      # 'Utility PV', 
       # 'IGCCAvgCF', 
+      # 'Biopower', 
+      'Land-Based Wind',
+      'Utility PV', 
       # 'LWR',
       # 'NuclearSMR', 
       # '10Hr Battery Storage', 
       # '8Hr Battery Storage',
       # '6Hr Battery Storage', 
-      # '4Hr Battery Storage', 
+      '4Hr Battery Storage', 
       # '2Hr Battery Storage'
       ]
 
diff --git a/dag.png b/dag.png
index 15b7047..89892c2 100644
Binary files a/dag.png and b/dag.png differ
diff --git a/functions/nrel_data_api.py b/functions/nrel_data_api.py
index f1aab75..9f6467e 100644
--- a/functions/nrel_data_api.py
+++ b/functions/nrel_data_api.py
@@ -21,7 +21,7 @@
               'year':2019,
               'leap_day':'false',
               'selector':'POINT',
-              'utc':'false',
+              'utc':'true',
               'interval':'60',
               'attr_list':['ghi']}
 
diff --git a/scripts/add_electricity.py b/scripts/add_electricity.py
index a60ecf2..c4fb0a2 100644
--- a/scripts/add_electricity.py
+++ b/scripts/add_electricity.py
@@ -6,19 +6,30 @@
 from tqdm import tqdm
 import pypsa
 
-model_year = snakemake.config['model_year']
+version = 22
+model_years = np.array(snakemake.config['model_years']).astype('int')
+resolution = int(snakemake.config['time_res'])
+scale = snakemake.config['geo_res']
+idx_opts = {"rto":"balancing_authority_code",
+            "county":"county"}
+growth_rates = snakemake.config['growth_rates']
+
+BUILD_YEAR = 2025  # a universal build year place holder
+
 
 def annuity(r, n):
     return r / (1-1/(1+r)**n)
 
 
 def create_snapshots():
-    snapshots = pd.date_range(start=f"{model_year}-01-01", 
-                            end=f"{model_year+1}-01-01", 
-                            freq='1h', 
-                            inclusive='left')
+    timestamps = pd.DatetimeIndex([])
+    for year in model_years:
+        period = pd.date_range(start=f"{year}-01-01",
+                                        freq=f"{resolution}h",
+                                        periods=8760/resolution)
+        timestamps = timestamps.append(period)
     
-    return snapshots
+    return timestamps
 
 
 def load_costs():
@@ -48,22 +59,64 @@ def load_costs():
 
 def load_existing_generators():
     generators = pd.read_csv(snakemake.input.generators, 
-                             index_col='balancing_authority_code')
+                             index_col=idx_opts[scale])
 
     return generators
+
+def load_build_years():
+    build_years = pd.read_csv(snakemake.input.build_years,
+                              index_col=idx_opts[scale])
     
+    return build_years
+    
+    
+def linear_growth(init_value, start_year, growth_rate, end_year=2050):
+    def model(x, init_val, start, rate):
+        return rate * init_val * (x - start) + init_val
+    years = np.arange(start_year, end_year, 1).astype('int')
+    growth_data = model(years, init_value, start_year, growth_rate)
+
+    growth_df = pd.DataFrame({'demand':growth_data})
+    growth_df.index = pd.date_range(start=str(start_year), 
+                                    periods=(end_year-start_year),
+                                    freq='YE')
 
-# attach components
+
+    return growth_df
+     
+
+# attach components 
 def attach_load(n):
-    load = pd.read_csv(snakemake.input.load, parse_dates=True, index_col="Interval End")
+    initial_demand = float(snakemake.config['total_demand'])
     
-    load = load.loc[str(snakemake.config['load_year'])]
+    if len(model_years) == 1:
+        demand = [initial_demand]
+    else:
+        rate = float(snakemake.config['load_growth'])
+        start_year = model_years[0]
+        end_year = model_years[-1]
+        N_years = end_year-start_year
+        demand = linear_growth(initial_demand, 
+                               model_years[0], 
+                               rate
+                               )
+        demand = demand.loc[demand.index.year.isin(model_years)].values.flatten()
+    
+    
+    load = pd.read_csv(snakemake.input.load, parse_dates=True, index_col="Interval End")
     
-    snapshots = create_snapshots()
+    start = int(snakemake.config['load_year'])
+    n_model_years = len(snakemake.config['model_years'])
+    end = start + n_model_years - 1
+    load = load.loc[str(start):str(end)][:int(n_model_years*8760)]
+    # normalize the load data
     
-    load.set_index(snapshots, inplace=True)
+    for d,y in zip(demand, load.index.year.unique()):
+        load.loc[str(y)] = load.loc[str(y)].div(load.loc[str(y)].sum(axis=0).sum(),axis=1)*d
     
     load = load.resample(f"{snakemake.config['time_res']}h").mean()
+    load = load[:len(n.snapshots)]
+    load.set_index(n.snapshots, inplace=True)
     
     print('Adding loads to model')
     for bus in tqdm(n.buses.index):
@@ -76,79 +129,194 @@ def attach_load(n):
     return
 
 
-def attach_generators(n, costs, generators):
-    costs = costs.xs((slice(None), slice(None), model_year))
-
-    available_carriers = costs.index.get_level_values('technology_alias').unique().to_list()
-    
-    snapshots = create_snapshots()
-    wind_profile = pd.read_csv(snakemake.input.wind_profile, parse_dates=True, index_col=0)
-    wind_profile.set_index(snapshots, inplace=True)
-    wind_profile = wind_profile.resample(f"{snakemake.config['time_res']}h").mean()
+def load_emissions():
+    df = pd.read_csv(snakemake.input.emissions, index_col=0)
     
-    solar_profile = pd.read_csv(snakemake.input.solar_profile, parse_dates=True, index_col=0)
-    solar_profile.set_index(snapshots, inplace=True)
-    solar_profile = solar_profile.resample(f"{snakemake.config['time_res']}h").mean()
+    return df
+
+
+def add_carriers(n, costs, emissions):
+    available_carriers = list(costs['technology_alias'].unique())
     
-    # add carriers
     for carrier in available_carriers:
+        # emissions data
+        try:
+            co2_emissions = emissions.at[carrier,'tCO2 per MWh']
+        except KeyError:
+            co2_emissions = 0.00
+            
         n.add(class_name="Carrier",
               name=carrier, 
-              color=snakemake.config['carrier_colors'][carrier])
-        
-    # flatten index
-    costs = costs.reset_index()
-    costs.loc[costs['technology_alias']=='Solar', 'techdetail'] = 'Utility PV'
+              co2_emissions=co2_emissions,
+              color=snakemake.config['carrier_colors'][carrier],
+              max_growth=float(growth_rates[carrier]))
+    return
+
+
+def load_re_profile(n, carrier='Solar'):
+    file_name = {'Solar':snakemake.input.solar_profile,
+                 'Wind':snakemake.input.wind_profile}  
+    re_profile = pd.read_csv(file_name[carrier], parse_dates=True, index_col=0)
+    re_profile = re_profile.resample(f"{snakemake.config['time_res']}h").mean().dropna(axis=0)
+    re_profile.set_index(n.snapshots, inplace=True) 
+    
+    return re_profile
+
+
+
+def attach_renewables(n, costs, generators=None, build_years=None, model_year=None):
+    carriers = ['Wind','Solar']
+    
     for bus in n.buses.index:
         for item in costs.itertuples():
             tech = item.techdetail
             carrier = item.technology_alias
-            capital_cost = item.capital_cost
-            marginal_cost = item.marginal_cost
-            lifetime = item.lifetime
-            
-            if carrier == 'Batteries':
-                class_name = "StorageUnit"
-                max_hours = float(tech.split(' ')[0].strip('Hr'))
-                cyclic_state_of_charge=True
-            else: 
-                class_name = "Generator"
-                max_hours = 0.0
-                cyclic_state_of_charge=False
-            
-            # existing capacity
-            if tech in generators.loc[bus].index:
-                p_nom = generators.at[bus, tech]
-            else:
-                p_nom = 0.0
-                
-            # renewables
-            if carrier == 'Wind':
-                p_max_pu = p_min_pu = wind_profile[bus]
-            elif carrier == 'Solar':
-                p_max_pu = p_min_pu = solar_profile[bus]
-            else:
-                p_max_pu = 1
-                p_min_pu = 0
-                
-            extendable = tech in snakemake.config['extendable_techs']
-
-            n.add(class_name=class_name,
-                  name=f"{bus} {tech}",
-                  bus=bus,
-                  p_nom=p_nom,
-                  p_nom_min=p_nom,
-                  p_min_pu=p_min_pu,
-                  p_max_pu=p_max_pu,
-                  p_nom_extendable=extendable,
-                  carrier=carrier,
-                  capital_cost=capital_cost,
-                  marginal_cost=marginal_cost,
-                  lifetime=lifetime,
-                  max_hours=max_hours,
-                  cyclic_state_of_charge=cyclic_state_of_charge)
+            if carrier in carriers:
+                re_profile = load_re_profile(n, carrier=carrier)
+                # existing capacity
+                if isinstance(generators, pd.DataFrame) and isinstance(build_years, pd.DataFrame):
+                    try:
+                        p_nom = generators.at[bus, tech]
+                        build_year = build_years.at[bus, tech]
+                    except:  # the technology is not included in the list of existing generators
+                        continue
+                    name = f"{bus} {tech} EXIST"
+                    extendable = False
+                elif model_year:
+                    build_year = model_year
+                    p_nom = 0.0
+                    name = f"{bus} {tech} {model_year}"
+                    extendable = tech in snakemake.config['extendable_techs']
+                    build_year = model_year
+                    
+                p_max_pu = re_profile[bus]
+                    
+
+                n.add(class_name="Generator",
+                    name=name,
+                    bus=bus,
+                    p_nom=p_nom,
+                    p_nom_min=p_nom,
+                    p_max_pu=p_max_pu,
+                    p_nom_extendable=extendable,
+                    carrier=carrier,
+                    capital_cost=item.capital_cost,
+                    marginal_cost=item.marginal_cost,
+                    lifetime=item.lifetime,
+                    build_year=build_year)
+    return 
+
+
+def attach_generators(n, costs, generators=None, build_years=None, model_year=None):
+    carriers = ['Natural Gas','Biomass','Coal','Nuclear','Petroleum']
+    for bus in n.buses.index:
+        for item in costs.itertuples():
+            tech = item.techdetail
+            carrier = item.technology_alias
+            if carrier in carriers:
+                # existing capacity
+                if isinstance(generators, pd.DataFrame) and isinstance(build_years, pd.DataFrame):
+                    try:
+                        p_nom = generators.at[bus, tech]
+                        build_year = build_years.at[bus, tech]
+                    except:  # the technology is not included in the list of existing generators
+                        continue
+                    name = f"{bus} {tech} EXIST"
+                    extendable = False
+                elif model_year:
+                    build_year = model_year
+                    p_nom = 0.0
+                    name = f"{bus} {tech} {model_year}"
+                    extendable = tech in snakemake.config['extendable_techs']
+                    build_year = model_year
+                    
+                # ramp limits
+                if tech in ['LWR']:
+                    ramp_limit_up = 0.01*resolution
+                    ramp_limit_down = 0.01*resolution
+                    p_nom_min = p_nom
+                elif tech in ['IGCCAvgCF', 'Biopower']:
+                    ramp_limit_up = 0.1*resolution
+                    ramp_limit_down = 0.1*resolution
+                    p_nom_min = 0
+                elif tech in ['CCAvgCF','CTAvgCF']:
+                    ramp_limit_up = min(0.6*resolution,1.0)
+                    ramp_limit_down = min(0.6*resolution,1.0)
+                    p_nom_min = 0
+                elif tech in ['Petroleum']:
+                    ramp_limit_up = 1.0
+                    ramp_limit_down = 1.0
+                    p_nom_min = 0
+                else:   
+                    ramp_limit_up = 1.0
+                    ramp_limit_down = 1.0
+                    p_nom_min = p_nom
+
+                # minimum/maximum power output
+                if tech == 'LWR':
+                    p_min_pu = 0.95
+                    p_max_pu = 1.0  
+                else:
+                    p_max_pu = 1
+                    p_min_pu = 0
+
+                n.add(class_name="Generator",
+                    name=name,
+                    bus=bus,
+                    p_nom=p_nom,
+                    p_nom_min=p_nom_min,
+                    p_nom_extendable=extendable,
+                    p_min_pu = p_min_pu,
+                    p_max_pu = p_max_pu,
+                    carrier=carrier,
+                    capital_cost=item.capital_cost,
+                    marginal_cost=item.marginal_cost,
+                    lifetime=item.lifetime,
+                    ramp_limit_down = ramp_limit_down,
+                    ramp_limit_up = ramp_limit_up,
+                    build_year=build_year,
+                    )
     return
-                
+              
+              
+def attach_storage(n, costs, generators=None, build_years=None, model_year=None):
+    carriers = ["Batteries"]
+    for bus in n.buses.index:
+        for item in costs.itertuples():
+            tech = item.techdetail
+            carrier = item.technology_alias
+            if carrier in carriers:
+                # existing capacity
+                if isinstance(generators, pd.DataFrame) and isinstance(build_years, pd.DataFrame):
+                    try:
+                        p_nom = generators.at[bus, tech]
+                        build_year = build_years.at[bus, tech]
+                    except:  # the technology is not included in the list of existing generators
+                        continue
+                    name = f"{bus} {tech} EXIST"
+                    extendable = False
+                elif model_year:
+                    build_year = model_year
+                    p_nom = 0.0
+                    name = f"{bus} {tech} {model_year}"
+                    extendable = tech in snakemake.config['extendable_techs']
+                    build_year = model_year
+
+                n.add(class_name="StorageUnit",
+                    name=name,
+                    bus=bus,
+                    p_nom=p_nom,
+                    p_nom_min=p_nom,
+                    p_nom_extendable=extendable,
+                    carrier=carrier,
+                    capital_cost=item.capital_cost,
+                    marginal_cost=item.marginal_cost,
+                    lifetime=item.lifetime,
+                    max_hours=float(tech.split(' ')[0].strip('Hr')),
+                    cyclic_state_of_charge=False,
+                    build_year=build_year)
+        
+    return  
 
 
 if __name__ == "__main__":
@@ -157,11 +325,64 @@ def attach_generators(n, costs, generators):
     costs = load_costs()
     
     costs.to_csv("data/final_costs.csv")
+    current_costs = costs.xs((slice(None), slice(None), 2020))
+    current_costs = current_costs.reset_index()
+    current_costs.loc[current_costs['technology_alias']=='Solar', 'techdetail'] = 'Utility PV'
     
     generators = load_existing_generators()
+    build_years = load_build_years()
+    emissions = load_emissions()
     
     attach_load(n)
-    attach_generators(n, costs=costs, generators=generators)    
+    add_carriers(n, 
+                 costs=current_costs,
+                 emissions=emissions)
+    
+    # add existing technology
+    attach_renewables(n,
+                      costs=current_costs,
+                      generators=generators,
+                      build_years=build_years
+                      )
+    attach_generators(n, 
+                    costs=current_costs, 
+                    generators=generators,
+                    build_years=build_years
+                    )   
+    attach_storage(n,
+                   costs=current_costs,
+                   generators=generators,
+                   build_years=build_years
+                   )
+    
+    # add new technology
+    for year in model_years:
+        # current_costs = costs.xs((slice(None), slice(None), year))
+        # current_costs = current_costs.reset_index()
+        # current_costs.loc[current_costs['technology_alias']=='Solar', 'techdetail'] = 'Utility PV'
+        attach_renewables(n,
+                        costs=current_costs,
+                        model_year=year
+                        )
+        attach_generators(n, 
+                        costs=current_costs, 
+                        model_year=year
+                        )   
+        attach_storage(n,
+                    costs=current_costs,
+                    model_year=year
+                    )
+    
+    # add co2 constraint    
+    for y in model_years:
+        emissions_dict = snakemake.config['co2_limits']
+        
+        n.add(class_name="GlobalConstraint",
+              name=f"CO2 Limit {y}",
+              carrier_attribute='co2_emissions',
+              sense="<=",
+              investment_period=y,
+              constant=float(emissions_dict[y])*1e6) 
     
     n.export_to_netcdf(snakemake.output.elec_network)
     
diff --git a/scripts/build_base_network.py b/scripts/build_base_network.py
index 95ff7b8..f1a7112 100644
--- a/scripts/build_base_network.py
+++ b/scripts/build_base_network.py
@@ -4,6 +4,27 @@
 import pypsa
 
 
+model_years = np.array(snakemake.config['model_years']).astype('int')
+resolution = int(snakemake.config['time_res'])
+
+def add_snapshots(network):
+
+    n_hours = 8760
+    
+    snapshots = pd.DatetimeIndex([])
+    for year in model_years:
+        period = pd.date_range(start=f"{year}-01-01",
+                                        freq=f"{resolution}h",
+                                        periods=n_hours / resolution)
+        snapshots = snapshots.append(period)
+    
+    network.snapshots = pd.MultiIndex.from_arrays([snapshots.year, snapshots])
+    network.investment_periods = model_years
+    
+    network.snapshot_weightings.loc[:,:] = resolution
+    
+    return
+
 def base_network():
     n = pypsa.Network()
     n.name = 'PyPSA-Illinois'
@@ -14,14 +35,7 @@ def base_network():
     line_config = snakemake.config['lines']
     v_nom = line_config['v_nom']
 
-    resolution = int(snakemake.config['time_res'])
-    model_year = snakemake.config['model_year']
-    n.set_snapshots(pd.date_range(start=f"{model_year}-01-01", 
-                                  end=f"{model_year+1}-01-01", 
-                                  inclusive='left',
-                                  freq=f"{resolution}h"))
     
-    n.snapshot_weightings.loc[:,:] = resolution
     n.import_components_from_dataframe(buses, 'Bus')
     n.import_components_from_dataframe(lines, 'Line')
     
@@ -29,4 +43,22 @@ def base_network():
 
 if __name__ == "__main__":
     n = base_network()
+    add_snapshots(n)
+    
+    year_diff = np.diff(model_years)
+    if len(year_diff) == 0:
+        pad = 1
+    else:
+        pad = year_diff[-1]
+    
+    n.investment_period_weightings["years"] = list(year_diff) + [pad]
+
+    r = 0.01
+    T = 0
+    for period, nyears in n.investment_period_weightings.years.items():
+        discounts = [(1 / (1 + r) ** t) for t in range(T, T + nyears)]
+        n.investment_period_weightings.at[period, "objective"] = sum(discounts)
+        T += nyears
+    n.investment_period_weightings
+    
     n.export_to_netcdf(snakemake.output.network)
\ No newline at end of file
diff --git a/scripts/plot_results.py b/scripts/plot_results.py
new file mode 100644
index 0000000..7a634dc
--- /dev/null
+++ b/scripts/plot_results.py
@@ -0,0 +1,149 @@
+import pypsa
+import pandas as pd
+import matplotlib.pyplot as plt
+import numpy as np
+
+def plot_dispatch(n, year=2025, month=7):
+    
+    time = (year, f'{year}-0{month}')
+    p_by_carrier = n.generators_t.p.groupby(n.generators.carrier, axis=1).sum().div(1e3)
+    
+    p_by_carrier = p_by_carrier[['Nuclear', 
+                                 'Coal', 
+                                 'Natural Gas', 
+                                 'Biomass', 
+                                 'Petroleum', 
+                                 'Solar',
+                                 'Wind']]
+
+    if not n.storage_units.empty:
+        sto = n.storage_units_t.p.T.groupby(n.storage_units.carrier).sum().T.div(1e3)
+        p_by_carrier = pd.concat([p_by_carrier, sto], axis=1)
+
+    fig, ax = plt.subplots(figsize=(12, 6))
+
+    color = p_by_carrier.columns.map(n.carriers.color)
+
+    p_by_carrier.where(p_by_carrier > 0).loc[time].plot.area(
+        ax=ax,
+        linewidth=0,
+        color=color,
+    )
+
+    charge = p_by_carrier.where(p_by_carrier < 0).dropna(how="all", axis=1).loc[time]
+
+    if not charge.empty:
+        charge.plot.area(
+            ax=ax,
+            linewidth=0,
+            color=charge.columns.map(n.carriers.color),
+        )
+
+    n.loads_t.p_set.sum(axis=1).loc[time].div(1e3).plot(ax=ax, c="k")
+
+    ax.legend(loc=(1.05, 0))
+    ax.set_ylabel("GW")
+    ax.set_ylim(-n.storage_units_t.p_store.max().max()/1e3-2.5, n.loads_t.p_set.sum(axis=1).max()/1e3+2.5)
+    plt.tight_layout()
+    
+    return fig, ax
+
+
+def plot_capacity(n):
+    
+    fig, ax = plt.subplots(figsize=(12,8))
+    
+    df = pd.concat([n.generators[['p_nom', 'p_nom_opt']], 
+                    n.storage_units[['p_nom','p_nom_opt']]])\
+                        .replace(0, np.nan)\
+                            .dropna(axis=0, how='all')
+                            
+    df.plot.bar(ax=ax)
+    plt.tight_layout()
+    
+    return fig,ax
+    
+    
+def plot_emissions(n, time_res):
+    
+    emissions_data = n.carriers\
+                        .reset_index()\
+                            .sort_values(by='Carrier')\
+                                .set_index('Carrier')[['co2_emissions']]\
+                                    .drop('Batteries')
+    
+    emissions = n.generators_t.p * n.generators.carrier.map(n.carriers.co2_emissions)
+    total_emissions = n.snapshot_weightings.generators @ emissions.sum(axis=1).div(1e6)
+    
+    p_by_carrier = n.generators_t.p.groupby(n.generators.carrier, axis=1).sum()
+    
+    p_by_carrier_year = (p_by_carrier.groupby(p_by_carrier.index.get_level_values('period')).sum())
+    
+    annual_emissions = (p_by_carrier_year * emissions_data['co2_emissions']).sum(axis=1).to_frame()
+    
+    annual_emissions = annual_emissions * time_res / 1e6
+    
+    annual_emissions.columns = ['Mtonnes CO2/year']
+
+    fig, ax = plt.subplots(figsize=(12,8))
+    
+    annual_emissions.plot.bar(ax=ax)
+    
+    ax.set_ylabel("Mtonnes CO2/year")
+    
+    return fig, ax
+
+
+def plot_active_units(n):
+    
+    fig, ax = plt.subplots(figsize=(12,8))
+    c = "StorageUnit"
+    df = pd.concat(
+        {
+            period: n.get_active_assets(c, period) * n.df(c).p_nom_opt
+            for period in n.investment_periods
+        },
+        axis=1,
+    )
+    df = df.groupby(n.storage_units.carrier).sum()
+
+    c = "Generator"
+    df2 = pd.concat(
+        {
+            period: n.get_active_assets(c, period) * n.df(c).p_nom_opt
+            for period in n.investment_periods
+        },
+        axis=1,
+    )
+    df2 = df2.groupby(n.generators.carrier).sum()
+    df = pd.concat([df, df2])
+    df.T.plot.bar(
+        ax=ax,
+        stacked=True,
+        edgecolor="white",
+        width=1,
+        ylabel="Capacity",
+        xlabel="Investment Period",
+        rot=0,
+        figsize=(10, 5),
+    )
+    plt.tight_layout()
+
+    return fig, ax
+
+if __name__ == "__main__":
+    
+    n = pypsa.Network(snakemake.input.solved_network)
+    
+    fig, ax = plot_dispatch(n, year=int(snakemake.config['plot_year']))
+    
+    plt.savefig(snakemake.output.dispatch_figure)
+    
+    fig, ax = plot_capacity(n)
+    plt.savefig(snakemake.output.capacity_figure)
+    
+    fig, ax = plot_emissions(n, float(snakemake.config['time_res']))
+    plt.savefig(snakemake.output.emissions_figure)
+    
+    fig, ax = plot_active_units(n)
+    plt.savefig(snakemake.output.active_units_figure)
\ No newline at end of file
diff --git a/scripts/retrieve_costs.py b/scripts/retrieve_costs.py
index d5b1d21..f4a93e1 100644
--- a/scripts/retrieve_costs.py
+++ b/scripts/retrieve_costs.py
@@ -1,6 +1,21 @@
 import pandas as pd
 from nrelpy.atb import ATBe
 
+n_illinois_reactors = 11
+total_lwr_capacity = 12415.1
+average_reactor_size = total_lwr_capacity/n_illinois_reactors  # MW
+nuclear_params = snakemake.config['nuclear_params']
+capacity_factor = float(nuclear_params['capacity_factor'])
+capital_cost = float(nuclear_params['capital_cost'])  # $/MWh
+nuclear_fuel = float(nuclear_params['fuel'])  # $/MWh
+fixed_om = float(nuclear_params['fixed_om'])  # $/MWh
+
+# convert capital and fom to $/MW/yr
+operating_hours = 8760*capacity_factor
+nuclear_cap_cost = capital_cost*operating_hours/1e3 # $/kW/yr
+nuclear_fixed_om = fixed_om*operating_hours/1e3  # $/kW/yr
+
+
 
 if __name__ == "__main__":
     atb_params = snakemake.config['atb_params']
@@ -85,11 +100,13 @@
     # add petroleum
     df_t = df_pivot.T
     
-    # these costs are from 2021 and should be updated to reflect inflation.
-    df_t['Petroleum','Petroleum',snakemake.config['model_year']] = [1158, 
-                                   27.94, 
-                                   1.78, 
-                                   petroleum_price*petroleum_heatrate]   # from here https://www.eia.gov/electricity/generatorcosts/
+    # [capital cost, fixed om cost, variable om cost, fuel cost]
+    for year in range(2020, 2051, 1):
+        # update costs for existing nuclear, these costs are from 2022
+        df_t['Nuclear','LWR',year] = [nuclear_cap_cost, nuclear_fixed_om, 0, nuclear_fuel]
+        
+        # these costs are from 2021 and should be updated to reflect inflation.
+        df_t['Petroleum','Petroleum',year] = [1158, 27.94, 1.78, petroleum_price*petroleum_heatrate]   # from here https://www.eia.gov/electricity/generatorcosts/
     df_pivot = df_t.T
     
     df_pivot.fillna(0., inplace=True)
diff --git a/scripts/retrieve_emissions.py b/scripts/retrieve_emissions.py
new file mode 100644
index 0000000..ba4da42
--- /dev/null
+++ b/scripts/retrieve_emissions.py
@@ -0,0 +1,18 @@
+import pandas as pd
+from unyt import g, kg, tonne, pound, kWh, MWh, GWh
+
+
+if __name__ == "__main__":
+
+    dfs = pd.read_html("https://www.eia.gov/tools/faqs/faq.php?id=74&t=11")
+    df = dfs[0].head(3)
+    df = df.replace({'Natural gas':'Natural Gas'})
+    df.columns = df.columns.droplevel()
+    df = df.set_index(df.iloc[:,0])
+    df = df[['pounds per kWh']]
+
+    df = df.iloc[:,0].astype('float').apply(lambda x: (x*(pound/kWh)).to(tonne/MWh)).to_frame()
+
+    df.index.name = ''
+    df.columns = ['tCO2 per MWh']
+    df.to_csv(snakemake.output.emissions)
\ No newline at end of file
diff --git a/scripts/retrieve_generators.py b/scripts/retrieve_generators.py
index be47b5f..a663854 100644
--- a/scripts/retrieve_generators.py
+++ b/scripts/retrieve_generators.py
@@ -126,10 +126,20 @@
     
     gdf.to_csv(snakemake.output.generators)
     
+    gdf['build_year'] = pd.to_datetime(gdf['operating-year-month']).dt.year
+    
     scale = snakemake.config['geo_res']
     idx_opts = {"rto":"balancing_authority_code",
                 "county":"county"}
     
+        
+    gen_build_year = gdf.pivot_table(index=idx_opts[scale],
+                                    columns='technology',
+                                    values='build_year',
+                                    aggfunc='mean').dropna(axis=1).astype('int')
+
+    gen_build_year.to_csv(snakemake.output.build_year)
+    
     gen_agg = gdf.pivot_table(index=idx_opts[scale],
                                 columns='technology',
                                 values='nameplate-capacity-mw',
diff --git a/scripts/retrieve_load.py b/scripts/retrieve_load.py
index 7a47997..d114645 100644
--- a/scripts/retrieve_load.py
+++ b/scripts/retrieve_load.py
@@ -30,6 +30,7 @@
                 verbose=True)
 
     try:
+        # demand['Interval End'] = demand['Interval End'].dt.tz_convert('US/Central')
         demand['Interval End'] = demand['Interval End'].dt.tz_localize(None)
     except:
         pass
@@ -40,7 +41,18 @@
                    columns='Subregion',
                    values='MW')
     
-    demand_pivot[demand_pivot>60e3] = np.nan # this should be replaced with a function to remove outliers. e.g., 2 SD away from mean.
+        
+    # calculate outliers with IQR
+    # threshold = 20
+    # q1 = demand_pivot.quantile(0.05)
+    # q3 = demand_pivot.quantile(0.95)
+    # IQR = q3-q1
+    # upper = (q3+threshold*IQR)
+    # lower = (q1-threshold*IQR)
+    # demand_pivot = demand_pivot[~((demand_pivot<lower) | (demand_pivot>upper))]
+    
+    
+    demand_pivot[demand_pivot>float(snakemake.config['load_filter'])] = np.nan # this should be replaced with a function to remove outliers
     demand_pivot = demand_pivot.interpolate("linear")
     
     demand_pivot.rename(columns=dict(zip(snakemake.config['rto_subba'],
diff --git a/scripts/retrieve_renewables.py b/scripts/retrieve_renewables.py
index 1e92c7c..a1043ce 100644
--- a/scripts/retrieve_renewables.py
+++ b/scripts/retrieve_renewables.py
@@ -3,12 +3,13 @@
 import pandas as pd
 import geopandas as gpd
 from tqdm import tqdm
-from unyt import m, s, MW, W, kg
+from unyt import m, s, MW, W, kg, g
 
 sys.path.append("functions")
 
 from nrel_data_api import parameters, make_csv_url
 
+model_years = np.array(snakemake.config['model_years']).astype('int')
 
 def handle_datetime(dataframe):
     """
@@ -21,10 +22,14 @@ def handle_datetime(dataframe):
         A pandas dataframe.
     """
     frame = dataframe.copy()
-    model_year = snakemake.config['model_year']
-    timestamps = pd.date_range(f"{model_year}-01-01",f"{model_year+1}-01-01",
-                               freq='1h', inclusive='left')
-    
+    timestamps = pd.DatetimeIndex([])
+    for year in model_years:
+        period = pd.date_range(start=f"{year}-01-01",
+                                        freq=f"1h",
+                                        periods=8760)
+        timestamps = timestamps.append(period)
+    
+    frame.index = timestamps
     try:
         frame.set_index(timestamps,inplace=True)
         frame.drop(columns=['Year','Month','Day','Hour','Minute'],inplace=True)
@@ -43,24 +48,32 @@ def retrieve_solar_timeseries(region):
     region : :class:`gpd.GeoDataFrame`
         A geopandas dataframe containing modeled bus regions.
     """
+   
     parameters['attr_list'] = ['ghi']
-    parameters['year'] = int(snakemake.config['solar_year'])
-    pbar = tqdm(region[['name','x','y']].values, position=0, leave=True)
-    frames = []
-    for n, i, j in pbar:
-        pbar.set_description(f"Processing {n}")
-        parameters['lon'] = i
-        parameters['lat'] = j
-        URL = make_csv_url(parameters=parameters, 
-                           kind='solar')
-        df = pd.read_csv(URL, skiprows=2)
-        df.rename(columns={'GHI':f"{n}"}, inplace=True)
-        df = handle_datetime(df)
-        frames.append(df)
-    
-    solar_df = pd.concat(frames, axis=1)
-    
-    return solar_df
+    outer_pbar = tqdm(snakemake.config['solar_years'], position=0, leave=True)
+    all_frames = []
+    for year in outer_pbar:
+        outer_pbar.set_description(f"Processing {year}")
+        parameters['year'] = int(year)
+        frames = []
+        inner_pbar = tqdm(region[['name','x','y']].values, position=1, leave=True)
+        for n, i, j in inner_pbar:
+            inner_pbar.set_description(f"Processing {n}")
+            parameters['lon'] = i
+            parameters['lat'] = j
+            URL = make_csv_url(parameters=parameters, 
+                            kind='solar')
+            df = pd.read_csv(URL, skiprows=2)[:8760]
+            df.rename(columns={'GHI':f"{n}"}, inplace=True)
+            frames.append(df)
+    
+            solar_df = pd.concat(frames, axis=1)
+    
+        all_frames.append(solar_df)
+    full_df = pd.concat(all_frames, axis=0)
+    full_df = handle_datetime(full_df)
+    
+    return full_df
 
 
 def retrieve_wind_timeseries(region):
@@ -74,23 +87,31 @@ def retrieve_wind_timeseries(region):
     """
     wind_attr = ['windspeed_80m']
     parameters['attr_list'] = wind_attr
-    parameters['year'] = int(snakemake.config['wind_year'])
-    pbar = tqdm(region[['name','x','y']].values, position=0, leave=True)
-    frames = []
-    for n, i, j in pbar:
-        pbar.set_description(f"Processing {n}")
-        parameters['lon'] = i
-        parameters['lat'] = j
-        URL = make_csv_url(parameters=parameters, 
-                           kind='wind')
-        df = pd.read_csv(URL, skiprows=1)
-        df.rename(columns={'wind speed at 80m (m/s)':f"{n}"}, inplace=True)
-        df = handle_datetime(df)
-        frames.append(df)
-    
-    wind_df = pd.concat(frames, axis=1)
-    
-    return wind_df
+    
+    outer_pbar = tqdm(snakemake.config['wind_years'], leave=True, position=0)
+    all_frames = []
+    for year in outer_pbar:
+        outer_pbar.set_description(f"Processing {year}")
+        
+        parameters['year'] = int(year)
+        pbar = tqdm(region[['name','x','y']].values, position=1, leave=True)
+        frames = []
+        for n, i, j in pbar:
+            pbar.set_description(f"Processing {n}")
+            parameters['lon'] = i
+            parameters['lat'] = j
+            URL = make_csv_url(parameters=parameters, 
+                            kind='wind')
+            df = pd.read_csv(URL, skiprows=1)[:8760]
+            df.rename(columns={'wind speed at 80m (m/s)':f"{n}"}, inplace=True)
+            frames.append(df)
+        
+        wind_df = pd.concat(frames, axis=1)
+        all_frames.append(wind_df)
+    full_df = pd.concat(all_frames, axis=0)
+    full_df = handle_datetime(full_df)
+    
+    return full_df
 
 
 def turbine_power(v):
diff --git a/scripts/solve_network.py b/scripts/solve_network.py
index f1fa404..7dd7c2e 100644
--- a/scripts/solve_network.py
+++ b/scripts/solve_network.py
@@ -3,58 +3,14 @@
 import pypsa
 import matplotlib.pyplot as plt
 
-def plot_dispatch(n, time="2025-07"):
-    p_by_carrier = n.generators_t.p.groupby(n.generators.carrier, axis=1).sum().div(1e3)
-    
-    p_by_carrier = p_by_carrier[['Nuclear', 
-                                 'Coal', 
-                                 'Natural Gas', 
-                                 'Biomass', 
-                                 'Petroleum', 
-                                 'Solar',
-                                 'Wind']]
 
-    if not n.storage_units.empty:
-        sto = n.storage_units_t.p.T.groupby(n.storage_units.carrier).sum().T.div(1e3)
-        p_by_carrier = pd.concat([p_by_carrier, sto], axis=1)
-
-    fig, ax = plt.subplots(figsize=(6, 3))
-
-    color = p_by_carrier.columns.map(n.carriers.color)
-
-    p_by_carrier.where(p_by_carrier > 0).loc[time].plot.area(
-        ax=ax,
-        linewidth=0,
-        color=color,
-    )
-
-    charge = p_by_carrier.where(p_by_carrier < 0).dropna(how="all", axis=1).loc[time]
-
-    if not charge.empty:
-        charge.plot.area(
-            ax=ax,
-            linewidth=0,
-            color=charge.columns.map(n.carriers.color),
-        )
-
-    n.loads_t.p_set.sum(axis=1).loc[time].div(1e3).plot(ax=ax, c="k")
+if __name__ == "__main__":
 
-    ax.legend(loc=(1.05, 0))
-    ax.set_ylabel("GW")
-    ax.set_ylim(n.generators_t.p.min().min(), n.loads_t.p_set.sum(axis=1).max()/1e3+2.5)
-    plt.tight_layout()
-    
-    return fig, ax
+    n = pypsa.Network(snakemake.input.elec_network)
 
+    try:
+        n.optimize(solver_name='cplex')
+    except:
+        n.optimize(solver_name='highs')
 
-if __name__ == "__main__":
-    
-    n = pypsa.Network(snakemake.input.elec_network)
-    
-    n.optimize(solver_name='highs')
-    
-    n.export_to_netcdf(snakemake.output.solved_network)
-    
-    fig, ax = plot_dispatch(n)
-    
-    plt.savefig(snakemake.output.dispatch_figure)
\ No newline at end of file
+    n.export_to_netcdf(snakemake.output.solved_network)
\ No newline at end of file