Negative Objective values observed in solved networks

[GbotemiB]

Hi, I am currently working on improving the Monte-Carlo feature in PyPSA-Earth.

Problem Objective

Making decisions is tough in Africa with fast-growing demand and big infrastructure changes. To deal with this uncertainty, I am using a Monte Carlo approach in the PyPSA-Earth model.

I am using Nigeria's Power System as a case study. I ran the simulation using the attached default config file with the Monte-Carlo option added, producing 961 networks. I observed that some of the solved networks had negative objective values.

My understanding of objective values translates to Annual cost which doesn't make sense to be negative.

# SPDX-FileCopyrightText:  PyPSA-Earth and PyPSA-Eur Authors
#
# SPDX-License-Identifier: CC0-1.0

version: 0.3.0
tutorial: false

logging:
  level: INFO
  format: "%(levelname)s:%(name)s:%(message)s"

countries: ["NG"]
# Can be replaced by country ["NG", "BJ"], continent ["Africa"] or user-specific region, see more at https://pypsa-earth.readthedocs.io/en/latest/configuration.html#top-level-configuration

enable:
  retrieve_databundle: false  #  Recommended 'true', for the first run. Otherwise data might be missing.
  retrieve_cost_data: false  # true: retrieves cost data from technology data and saves in resources/costs.csv, false: uses cost data in data/costs.csv
  download_osm_data: false  # If 'true', OpenStreetMap data will be downloaded for the above given countries
  build_natura_raster: false  # If True, then an exclusion raster will be build
  build_cutout: false
  # If "build_cutout" : true, then environmental data is extracted according to `snapshots` date range and `countries`
  # requires cds API key https://cds.climate.copernicus.eu/api-how-to
  # More information https://atlite.readthedocs.io/en/latest/introduction.html#datasets

custom_rules: []  # Default empty [] or link to custom rule file e.g. ["my_folder/my_rules.smk"] that add rules to Snakefile

run:
  name: "NG" # use this to keep track of runs with different settings
  shared_cutouts: true  # set to true to share the default cutout(s) across runs
                        # Note: value false requires build_cutout to be enabled

scenario:
  simpl: ['']
  ll: ['copt']
  clusters: [10]
  opts: [Co2L-24H]

snapshots:
  start: "2013-01-01"
  end: "2014-01-01"
blob:vscode-webview://0opju7i2n0lkb292tln19reu4jtuo9eg41lqjna9qh6plhob0bfe/abfb4ecf-09d3-4bfc-9e3a-515a135f56fe   inclusive: "left" # end is not inclusive

# definition of the Coordinate Reference Systems
crs:
  geo_crs: EPSG:4326  # general geographic projection, not used for metric measures. "EPSG:4326" is the standard used by OSM and google maps
  distance_crs: EPSG:3857  # projection for distance measurements only. Possible recommended values are "EPSG:3857" (used by OSM and Google Maps)
  area_crs: ESRI:54009  # projection for area measurements only. Possible recommended values are Global Mollweide "ESRI:54009"

# download_osm_data_nprocesses: 10  # (optional) number of threads used to download osm data

augmented_line_connection:
  add_to_snakefile: false  # If True, includes this rule to the workflow
  connectivity_upgrade: 2  # Min. lines connection per node, https://networkx.org/documentation/stable/reference/algorithms/generated/networkx.algorithms.connectivity.edge_augmentation.k_edge_augmentation.html#networkx.algorithms.connectivity.edge_augmentation.k_edge_augmentation
  new_line_type: ["HVAC"]  # Expanded lines can be either ["HVAC"] or ["HVDC"] or both ["HVAC", "HVDC"]
  min_expansion: 1  # [MW] New created line expands by float/int input
  min_DC_length: 600  # [km] Minimum line length of DC line

cluster_options:
  simplify_network:
    to_substations: false # network is simplified to nodes with positive or negative power injection (i.e. substations or offwind connections)
    algorithm: kmeans # choose from: [hac, kmeans]
    feature: solar+onwind-time # only for hac. choose from: [solar+onwind-time, solar+onwind-cap, solar-time, solar-cap, solar+offwind-cap] etc.
    exclude_carriers: []
    remove_stubs: true
    remove_stubs_across_borders: true
    p_threshold_drop_isolated: 20 # [MW] isolated buses are being discarded if bus mean power is below the specified threshold
    p_threshold_merge_isolated: 300 # [MW] isolated buses are being merged into a single isolated bus if bus mean power is below the specified threshold
  cluster_network:
    algorithm: kmeans
    feature: solar+onwind-time
    exclude_carriers: []
  alternative_clustering: false  # "False" use Voronoi shapes, "True" use GADM shapes
  distribute_cluster: ['load'] # Distributes cluster nodes per country according to ['load'],['pop'] or ['gdp']
  out_logging: true  # When "True", logging is printed to console
  aggregation_strategies:
    generators:  # use "min" for more conservative assumptions
      p_nom: sum
      p_nom_max: sum
      p_nom_min: sum
      p_min_pu: mean
      marginal_cost: mean
      committable: any
      ramp_limit_up: max
      ramp_limit_down: max
      efficiency: mean

build_shape_options:
  gadm_layer_id: 1  # GADM level area used for the gadm_shapes. Codes are country-dependent but roughly: 0: country, 1: region/county-like, 2: municipality-like
  update_file: false  # When true, all the input files are downloaded again and replace the existing files
  out_logging: true  # When true, logging is printed to console
  year: 2020  # reference year used to derive shapes, info on population and info on GDP
  nprocesses: 3  # number of processes to be used in build_shapes
  worldpop_method: "standard"  # "standard" pulls from web 1kmx1km raster, "api" pulls from API 100mx100m raster, false (not "false") no pop addition to shape which is useful when generating only cutout
  gdp_method: "standard"  # "standard" pulls from web 1x1km raster, false (not "false") no gdp addition to shape which useful when generating only cutout
  contended_flag: "set_by_country" # "set_by_country" assigns the contended areas to the countries according to the GADM database, "drop" drops these contended areas from the model

clean_osm_data_options:  # osm = OpenStreetMap
  names_by_shapes: true  # Set the country name based on the extended country shapes
  threshold_voltage: 51000  # [V] assets below that voltage threshold will not be used (cable, line, generator, etc.)
  tag_substation: "transmission"  # Filters only substations with 'transmission' tag, ('distribution' also available)
  add_line_endings: true  # When "True", then line endings are added to the dataset of the substations
  generator_name_method: OSM  # Methodology to specify the name to the generator. Options: OSM (name as by OSM dataset), closest_city (name by the closest city)
  use_custom_lines: "OSM_only"  # Use OSM (OSM_only), customized (custom_only), or both data sets (add_custom)
  path_custom_lines: false  # If exists, provide the specific absolute path of the custom file e.g. (...\data\custom_lines.geojson)
  use_custom_substations: "OSM_only"  # Use OSM (OSM_only), customized (custom_only), or both data sets (add_custom)
  path_custom_substations: false  # If exists, provide the specific absolute path of the custom file e.g. (...\data\custom_substations.geojson)
  use_custom_cables: "OSM_only"  # Use OSM (OSM_only), customized (custom_only), or both data sets (add_custom)
  path_custom_cables: false  # If exists, provide the specific absolute path of the custom file e.g. (...\data\custom_cables.geojson)

build_osm_network:  # Options of the build_osm_network script; osm = OpenStreetMap
  group_close_buses: true  # When "True", close buses are merged and guarantee the voltage matching among line endings
  group_tolerance_buses: 5000  # [m] (default 5000) Tolerance in meters of the close buses to merge
  split_overpassing_lines: true  # When True, lines overpassing buses are splitted and connected to the bueses
  overpassing_lines_tolerance: 1  # [m] (default 1) Tolerance to identify lines overpassing buses
  force_ac: false  # When true, it forces all components (lines and substation) to be AC-only. To be used if DC assets create problem.

base_network:
  min_voltage_substation_offshore: 51000  # [V] minimum voltage of the offshore substations
  min_voltage_rebase_voltage: 51000  # [V] minimum voltage in base network

load_options:
  ssp: "ssp2-2.6" # shared socio-economic pathway (GDP and population growth) scenario to consider
  weather_year: 2013  # Load scenarios available with different weather year (different renewable potentials)
  prediction_year: 2030  # Load scenarios available with different prediction year (GDP, population)
  scale: 1  # scales all load time-series, i.e. 2 = doubles load

electricity:
  base_voltage: 380.
  voltages: [220., 300., 380.]
  co2limit: 7.75e+7  # European default, 0.05 * 3.1e9*0.5, needs to be adjusted for Africa
  co2base: 1.487e+9  # European default, adjustment to Africa necessary
  agg_p_nom_limits: data/agg_p_nom_minmax.csv
  hvdc_as_lines: false  # should HVDC lines be modeled as `Line` or as `Link` component?
  automatic_emission: false
  automatic_emission_base_year: 1990 # 1990 is taken as default. Any year from 1970 to 2018 can be selected.

  operational_reserve: # like https://genxproject.github.io/GenX/dev/core/#Reserves
    activate: false
    epsilon_load: 0.02 # share of total load
    epsilon_vres: 0.02 # share of total renewable supply
    contingency: 0 # fixed capacity in MW

  max_hours:
    battery: 6
    H2: 168

  extendable_carriers:
    Generator: [solar, onwind, offwind-ac, offwind-dc, OCGT]
    StorageUnit: []  # battery, H2
    Store: [battery, H2]
    Link: []  # H2 pipeline

  powerplants_filter: (DateOut >= 2022 or DateOut != DateOut)
  custom_powerplants: false  #  "false" use only powerplantmatching (ppm) data, "merge" combines ppm and custom powerplants, "replace" use only custom powerplants

  conventional_carriers: [nuclear, oil, OCGT, CCGT, coal, lignite, geothermal, biomass]
  renewable_carriers: [solar, onwind, offwind-ac, offwind-dc, hydro]

  estimate_renewable_capacities:
    stats: "irena"  # False, = greenfield expansion, 'irena' uses IRENA stats to add expansion limits
    year: 2020  # Reference year, available years for IRENA stats are 2000 to 2020
    p_nom_min: 1  # any float, scales the minimum expansion acquired from stats, i.e. 110% of <years>'s capacities => p_nom_min: 1.1
    p_nom_max: false  # sets the expansion constraint, False to deactivate this option and use estimated renewable potentials determine by the workflow, float scales the p_nom_min factor accordingly
    technology_mapping:
      # Wind is the Fueltype in ppm.data.Capacity_stats, onwind, offwind-{ac,dc} the carrier in PyPSA-Earth
      Offshore: [offwind-ac, offwind-dc]
      Onshore: [onwind]
      PV: [solar]

lines:
  types:
    220.: "Al/St 240/40 2-bundle 220.0"
    300.: "Al/St 240/40 3-bundle 300.0"
    380.: "Al/St 240/40 4-bundle 380.0"
  dc_type: "HVDC XLPE 1000"
  s_max_pu: 0.7
  s_nom_max: .inf
  length_factor: 1.25
  under_construction: "zero" # 'zero': set capacity to zero, 'remove': remove, 'keep': with full capacity

links:
  p_max_pu: 1.0
  p_nom_max: .inf
  under_construction: "zero" # 'zero': set capacity to zero, 'remove': remove, 'keep': with full capacity

transformers:
  x: 0.1
  s_nom: 2000.
  type: ""

atlite:
  nprocesses: 4
  cutouts:
    # geographical bounds automatically determined from countries input
    cutout-2013-era5:
      module: era5
      dx: 0.3  # cutout resolution
      dy: 0.3  # cutout resolution
      # The cutout time is automatically set by the snapshot range. See `snapshot:` option above and 'build_cutout.py'.
      # time: ["2013-01-01", "2014-01-01"]  # to manually specify a different weather year (~70 years available)
      # The cutout spatial extent [x,y] is automatically set by country selection. See `countires:` option above and 'build_cutout.py'.
      # x: [-12., 35.]  # set cutout range manual, instead of automatic by boundaries of country
      # y: [33., 72]    # manual set cutout range


renewable:
  onwind:
    cutout: cutout-2013-era5
    resource:
      method: wind
      turbine: Vestas_V112_3MW
    capacity_per_sqkm: 3 # conservative, ScholzPhd Tab 4.3.1: 10MW/km^2
    # correction_factor: 0.93
    copernicus:
      # Scholz, Y. (2012). Renewable energy based electricity supply at low costs:
      #  development of the REMix model and application for Europe. ( p.42 / p.28)
      # CLC grid codes:
      # 11X/12X - Various forest types
      # 20  - Shrubs
      # 30  - Herbaceus vegetation
      # 40  - Cropland
      # 50  - Urban
      # 60  - Bare / Sparse vegetation
      # 80  - Permanent water bodies
      # 100 - Moss and lichen
      # 200 - Open sea
      grid_codes: [20, 30, 40, 60, 100, 111, 112, 113, 114, 115, 116, 121, 122, 123, 124, 125, 126]
      distance: 1000
      distance_grid_codes: [50]
    natura: true
    potential: simple # or conservative
    clip_p_max_pu: 1.e-2
    extendable: true
  offwind-ac:
    cutout: cutout-2013-era5
    resource:
      method: wind
      turbine: NREL_ReferenceTurbine_5MW_offshore
    capacity_per_sqkm: 2
    correction_factor: 0.8855
    # proxy for wake losses
    # from 10.1016/j.energy.2018.08.153
    # until done more rigorously in #153
    copernicus:
      grid_codes: [80, 200]
    natura: true
    max_depth: 50
    max_shore_distance: 30000
    potential: simple # or conservative
    clip_p_max_pu: 1.e-2
    extendable: true
  offwind-dc:
    cutout: cutout-2013-era5
    resource:
      method: wind
      turbine: NREL_ReferenceTurbine_5MW_offshore
    # ScholzPhd Tab 4.3.1: 10MW/km^2
    capacity_per_sqkm: 3
    correction_factor: 0.8855
    # proxy for wake losses
    # from 10.1016/j.energy.2018.08.153
    # until done more rigorously in #153
    copernicus:
      grid_codes: [80, 200]
    natura: true
    max_depth: 50
    min_shore_distance: 30000
    potential: simple # or conservative
    clip_p_max_pu: 1.e-2
    extendable: true
  solar:
    cutout: cutout-2013-era5
    resource:
      method: pv
      panel: CSi
      orientation: latitude_optimal # will lead into optimal design
        # slope: 0.  # slope: 0 represent a flat panel
        # azimuth: 180.  # azimuth: 180 south orientation
    capacity_per_sqkm: 4.6 # From 1.7 to 4.6 addresses issue #361
    # Determined by comparing uncorrected area-weighted full-load hours to those
    # published in Supplementary Data to
    # Pietzcker, Robert Carl, et al. "Using the sun to decarbonize the power
    # sector: The economic potential of photovoltaics and concentrating solar
    # power." Applied Energy 135 (2014): 704-720.
    correction_factor: 0.854337
    copernicus:
      grid_codes: [20, 30, 40, 50, 60, 90, 100]
    natura: true
    potential: simple # or conservative
    clip_p_max_pu: 1.e-2
    extendable: true
  hydro:
    cutout: cutout-2013-era5
    hydrobasins_level: 6
    resource:
      method: hydro
      hydrobasins: data/hydrobasins/hybas_world.shp
      flowspeed: 1.0  # m/s
      # weight_with_height: false
      # show_progress: true
    carriers: [ror, PHS, hydro]
    PHS_max_hours: 6
    hydro_max_hours: "energy_capacity_totals_by_country"  # not active
    hydro_max_hours_default: 6.0  # (optional, default 6) Default value of max_hours for hydro when NaN values are found
    clip_min_inflow: 1.0
    extendable: true
    normalization:
      method: hydro_capacities  # 'hydro_capacities' to rescale country hydro production by using hydro_capacities, 'eia' to rescale by eia data, false for no rescaling
      year: 2013  # (optional) year of statistics used to rescale the runoff time series. When not provided, the weather year of the snapshots is used
    multiplier: 1.1  # multiplier applied after the normalization of the hydro production; default 1.0

# TODO: Needs to be adjusted for Africa.
# Costs Configuration (Do not remove, needed for Sphynx documentation).
costs:
  year: 2030
  version: v0.5.0
  rooftop_share: 0.14  # based on the potentials, assuming  (0.1 kW/m2 and 10 m2/person)
  USD2013_to_EUR2013: 0.7532 # [EUR/USD] ECB: https://www.ecb.europa.eu/stats/exchange/eurofxref/html/eurofxref-graph-usd.en.html
  fill_values:
    FOM: 0
    VOM: 0
    efficiency: 1
    fuel: 0
    investment: 0
    lifetime: 25
    CO2 intensity: 0
    discount rate: 0.07
  marginal_cost: # EUR/MWh
    solar: 0.01
    onwind: 0.015
    offwind: 0.015
    hydro: 0.
    H2: 0.
    electrolysis: 0.
    fuel cell: 0.
    battery: 0.
    battery inverter: 0.
  emission_prices: # in currency per tonne emission, only used with the option Ep
    co2: 0.


monte_carlo:
    # Uniform: https://chaospy.readthedocs.io/en/master/api/chaospy.Uniform.html
    # Normal: https://chaospy.readthedocs.io/en/master/api/chaospy.Normal.html
    # LogNormal: https://chaospy.readthedocs.io/en/master/api/chaospy.LogNormal.html
    # Triange: https://chaospy.readthedocs.io/en/master/api/chaospy.Triangle.html
    # Beta: https://chaospy.readthedocs.io/en/master/api/chaospy.Beta.html
    # Gamma: https://chaospy.readthedocs.io/en/master/api/chaospy.Gamma.html
    # [mean, std] for Normal
    # [shape] for LogNormal
    # [lower_bound, upper_bound] for Uniform
    # [lower_bound, midpoint, upper_bound] for Triangle
    # [alpha, beta] for Beta
    # [shape, scale] for Gamma
    # TODO: Support inputs to simulate outages biggest lines "lines.p_nom_opt.max(): [-3000MW 0MW]""
    # TODO: Support inputs to simulate outages of biggest power plant "generators.p_nom.max()": [-1000MW 0MW]
    # ... user can add flexibly more features for the Monte-Carlo sampling
    # as key: add below the pypsa object for the monte_carlo sampling, "network" is only allowed for filtering!
    # as value: currently supported format [l_bound, u_bound] or empty [], represent multiplication factors for the object
  options:
    add_to_snakefile: true
    samples: 961 # number of optimizations. Note that number of samples when using scipy has to be the square of a prime number
    sampling_strategy: "chaospy"  # "pydoe2", "chaospy", "scipy", packages that are supported
    seed: 42 # set seedling for reproducibilty
  uncertainties:
    loads_t.p_set:
      type: beta
      args: [2, 0.5]
    generators_t.p_max_pu.loc[:, n.generators.carrier == "solar"]:
      type: beta
      args: [4, 1]
    generators_t.p_max_pu.loc[:, n.generators.carrier == "onwind"]:
      type: beta
      args: [4, 1]
    generators.capital_cost.loc[n.generators.carrier == "solar"]:
      type: normal
      args: [0, 1]
    generators.capital_cost.loc[n.generators.carrier == "onwind"]:
      type: normal
      args: [0, 1]
    generators.weight.loc[n.generators.carrier == "onwind"]:
      type: beta
      args: [2, 0.5]
    generators.weight.loc[n.generators.carrier == "solar"]:
      type: beta
      args: [2, 0.5]

solving:
  options:
    formulation: kirchhoff
    load_shedding: true
    noisy_costs: true
    min_iterations: 4
    max_iterations: 6
    clip_p_max_pu: 0.01
    skip_iterations: true
    track_iterations: false
    #nhours: 10
  solver:
    name: glpk
    # threads: 1
    # method: 2 # barrier (=ipm)
    # crossover: 0
    # BarConvTol: 1.e-5
    # FeasibilityTol: 1.e-6
    # AggFill: 0
    # PreDual: 0
    # GURO_PAR_BARDENSETHRESH: 200


plotting:
  map:
    figsize: [7, 7]
    boundaries: [-10.2, 29, 35, 72]
    p_nom:
      bus_size_factor: 5.e+4
      linewidth_factor: 3.e+3

  costs_max: 800
  costs_threshold: 1

  energy_max: 15000.
  energy_min: -10000.
  energy_threshold: 50.

  vre_techs: ["onwind", "offwind-ac", "offwind-dc", "solar", "ror"]
  conv_techs: ["OCGT", "CCGT", "nuclear", "coal", "oil"]
  storage_techs: ["hydro+PHS", "battery", "H2"]
  load_carriers: ["AC load"]
  AC_carriers: ["AC line", "AC transformer"]
  link_carriers: ["DC line", "Converter AC-DC"]
  tech_colors:
    "onwind": "#235ebc"
    "onshore wind": "#235ebc"
    "offwind": "#6895dd"
    "offwind-ac": "#6895dd"
    "offshore wind": "#6895dd"
    "offshore wind ac": "#6895dd"
    "offwind-dc": "#74c6f2"
    "offshore wind dc": "#74c6f2"
    "hydro": "#08ad97"
    "hydro+PHS": "#08ad97"
    "PHS": "#08ad97"
    "hydro reservoir": "#08ad97"
    "hydroelectricity": "#08ad97"
    "ror": "#4adbc8"
    "run of river": "#4adbc8"
    "solar": "#f9d002"
    "solar PV": "#f9d002"
    "solar thermal": "#ffef60"
    "biomass": "#0c6013"
    "solid biomass": "#06540d"
    "biogas": "#23932d"
    "waste": "#68896b"
    "geothermal": "#ba91b1"
    "OCGT": "#d35050"
    "gas": "#d35050"
    "natural gas": "#d35050"
    "CCGT": "#b20101"
    "nuclear": "#ff9000"
    "coal": "#707070"
    "lignite": "#9e5a01"
    "oil": "#262626"
    "H2": "#ea048a"
    "hydrogen storage": "#ea048a"
    "battery": "#b8ea04"
    "Electric load": "#f9d002"
    "electricity": "#f9d002"
    "lines": "#70af1d"
    "transmission lines": "#70af1d"
    "AC-AC": "#70af1d"
    "AC line": "#70af1d"
    "links": "#8a1caf"
    "HVDC links": "#8a1caf"
    "DC-DC": "#8a1caf"
    "DC link": "#8a1caf"
    "load": "#FF0000"
  nice_names:
    OCGT: "Open-Cycle Gas"
    CCGT: "Combined-Cycle Gas"
    offwind-ac: "Offshore Wind (AC)"
    offwind-dc: "Offshore Wind (DC)"
    onwind: "Onshore Wind"
    solar: "Solar"
    PHS: "Pumped Hydro Storage"
    hydro: "Reservoir & Dam"
    battery: "Battery Storage"
    H2: "Hydrogen Storage"
    lines: "Transmission Lines"
    ror: "Run of River"

Here is a plot of the objective values for all 961 solved networks

Here are the distribution plots for the parameters to which uncertainty is applied.

[ekatef]

According to the discussion today, the most likely reason for the negative objective is that the optimal value nom_opt of the parameter gets lower than minimal value nom_min. #956 provides the details for that. Thanks a lot @davide-f and @pz-max for the insights!

@GbotemiB thank you so much for sharing the results. feel free to check this explanation and propose a fix 😉