Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Power Portfolio Build-Out Example#

This notebook provides a simplified example of building a power generation portfolio using linopy. The goal is to demonstrate the process of setting up and solving a linear optimization problem where different power generation resources are selected and built to meet future energy demands.

[1]:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr

import linopy as lp

Problem Statement#

This notebook presents an example portfolio expansion model used to meet an energy need. This was originally developed for a hydropower utility. Hydropower utilities can be energy limited (i.e., monthly energy balance) rather than capacity limited (i.e., hourly generating capacity), therefore, this algorithm addresses a monthly energy need rather than an hourly capacity need. A different algorithm would be needed to meet capacity needs and would include battery storage, etc.

This notebook addresses the question: what is the least-cost way to meet an energy need in the future?

Generate Dummy Data#

The following cells only generate synthetic input data for the example. You can skim them and focus on the linopy model, which is defined further below. ### Energy Need

This analysis assumes that an energy need has been identified through modeling and simulation. For this example, we will assume the energy need is as follows:

[2]:
def generate_dummy_energy_needs(start_year=2024, years=10):
    # Generate a monthly date range
    dates = pd.date_range(start=f"{start_year}-01-01", periods=years * 12, freq="ME")

    # Create a sinusoidal pattern for energy needs with a downward trend over time
    months = np.arange(len(dates))

    # Base seasonal pattern using sine functions (with phase shift for winter, spring, etc.)
    seasonal_pattern = np.sin(
        2 * np.pi * (months % 12) / 12
    )  # Basic seasonal pattern (positive in spring/fall, negative in winter/summer)

    # Adjust the pattern so that winter months trend more negative over time
    winter_weight = 0.5 * (
        np.cos(2 * np.pi * (months % 12) / 12 - np.pi) + 1
    )  # Emphasize winter more

    # Add a downward trend over time
    trend = -0.01 * months  # A small negative trend each month

    # Combine seasonal pattern, winter weight, and trend
    energy_balance = seasonal_pattern * (1 - winter_weight) + trend
    energy_balance *= 10
    # Convert to pandas DataFrame
    df = pd.DataFrame(data={"monthly_energy_balance": energy_balance}, index=dates)

    return df
[3]:
# Generate the dummy dataset
dummy_energy_needs_df = generate_dummy_energy_needs()
dummy_energy_needs_df.plot()
plt.grid()
_images/energy-resource-build_5_0.svg

The energy need is the shortfall — the negative part of the monthly energy balance, expressed in average megawatts (aMW). It becomes an input to the model, so we store it as an xarray DataArray.

[4]:
dummy_energy_needs_df.index.name = "datetime"
dummy_energy_needs_df["energy_needs_aMW"] = (
    -dummy_energy_needs_df.monthly_energy_balance
).clip(lower=0)
energy_needs_xr = xr.DataArray(dummy_energy_needs_df.energy_needs_aMW)

Resource Generation Profiles#

Each resource has a normalized generation profile — its output per unit of installed capacity in each month. We make up dummy monthly profiles for wind and solar at three sites (A, B, C).

[5]:
def generate_resource_profiles(start_year=2024, years=10, noise_level=0.05):
    # Generate a monthly date range
    dates = pd.date_range(start=f"{start_year}-01-01", periods=years * 12, freq="ME")

    # Base profiles for wind and solar resources
    wind_winter_peak = np.array(
        [0.6, 0.55, 0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.3, 0.4, 0.5, 0.55]
    )  # Wind peaking in winter
    wind_summer_peak = np.array(
        [0.3, 0.25, 0.3, 0.4, 0.5, 0.55, 0.6, 0.55, 0.5, 0.45, 0.4, 0.35]
    )  # Wind peaking in summer
    solar_base = np.array(
        [0.1, 0.2, 0.4, 0.6, 0.7, 0.9, 0.95, 0.85, 0.6, 0.4, 0.2, 0.1]
    )  # Solar peaking in summer

    # Create dummy profiles with slight variations
    profiles = {}
    for i, location in enumerate(["A", "B", "C"]):
        # Seed a separate generator per location so the profiles stay reproducible
        wind_seed = hash(f"wind_{location}") % (2**32)
        solar_seed = hash(f"solar_{location}") % (2**32)

        wind_rng = np.random.default_rng(wind_seed)
        if i == 0:
            wind_profile = wind_winter_peak + noise_level * wind_rng.standard_normal(12)
        elif i == 1:
            wind_profile = wind_summer_peak + noise_level * wind_rng.standard_normal(12)
        else:
            # Create a mix or a balanced profile as the third option
            wind_profile = 0.5 * (
                wind_winter_peak + wind_summer_peak
            ) + noise_level * wind_rng.standard_normal(12)

        solar_rng = np.random.default_rng(solar_seed)
        solar_profile = solar_base + noise_level * solar_rng.standard_normal(12)

        # Separate generators for the final noise adder
        wind_noise_rng = np.random.default_rng(wind_seed + 1)
        profiles[f"wind_{location}"] = np.tile(
            wind_profile, years
        ) + noise_level * wind_noise_rng.standard_normal(len(dates))

        solar_noise_rng = np.random.default_rng(solar_seed + 1)
        profiles[f"solar_{location}"] = np.tile(
            solar_profile, years
        ) + noise_level * solar_noise_rng.standard_normal(len(dates))

    # Create a DataFrame with the profiles
    df_profiles = pd.DataFrame(data=profiles, index=dates)

    return df_profiles
[6]:
# Generate the dummy dataset
df_profiles = generate_resource_profiles()

# Plot the profiles for the first 12 months
df_profiles.groupby(df_profiles.index.month).mean().plot()
plt.xlabel("month")
plt.ylabel("normalized generation")
plt.grid()
_images/energy-resource-build_10_0.svg

As before, convert to xarray for linopy.

[7]:
df_profiles.index.name = "datetime"
df_profiles.columns.name = "resource"
resource_profiles_xr = xr.DataArray(df_profiles)

Resource Cost Profiles#

Each resource has a unit cost for building capacity. Here the costs decline slowly over time, reflecting technology learning — which gives the model a reason to weigh building sooner against building later.

[8]:
def generate_unit_cost_profiles(start_year=2024, years=10, decay_rate=0.005):
    # Generate a monthly date range
    dates = pd.date_range(start=f"{start_year}-01-01", periods=years * 12, freq="ME")

    # Initial costs for wind and solar resources
    initial_costs = {
        "wind_A": 1.0,
        "wind_B": 1.2,
        "wind_C": 1.1,
        "solar_A": 0.8,
        "solar_B": 0.85,
        "solar_C": 0.9,
    }

    # Create dummy profiles with slow exponential decay
    profiles = {}
    for resource, initial_cost in initial_costs.items():
        # Exponential decay function
        decay = initial_cost * np.exp(-decay_rate * np.arange(len(dates)))

        # Assign decay curve to profile
        profiles[resource] = decay

    # Create a DataFrame with the profiles
    df_profiles = pd.DataFrame(data=profiles, index=dates)

    return df_profiles


# Generate the dummy dataset
df_unit_cost_profiles = generate_unit_cost_profiles()

# Plot the unit cost profiles
df_unit_cost_profiles.plot()
plt.xlabel("dt")
plt.ylabel("Unit Cost")
plt.legend(title="resource", bbox_to_anchor=(1.05, 1), loc="upper left")
plt.grid()
_images/energy-resource-build_14_0.svg

Convert to xarray

[9]:
df_unit_cost_profiles.index.name = "datetime"
df_unit_cost_profiles.columns.name = "resource"
unit_cost_profiles_xr = xr.DataArray(df_unit_cost_profiles)

Optimized Energy Portfolio Model#

With the energy need, generation profiles, and costs in hand, we can write the optimization model. The question it answers is: how much of each resource should we build, and when, to meet the energy need at least cost?

The model has three ingredients, added to an empty linopy.Model across the cells below:

  • a decision variable, build_out, holding the installed capacity of each resource in each month,

  • constraints that keep the build-out physically sensible and force it to cover the energy need, and

  • an objective that minimizes the total build cost.

Because the variable is indexed over both datetime and resource, a single line of linopy expands to a whole array of variables or constraints.

[10]:
m = lp.Model()

The build-out variable#

build_out is the model’s single decision variable, indexed over datetime and resource: build_out[t, r] is the cumulative installed capacity of resource r in month t, bounded below by zero.

Capacity can be added over time but never retired, so the build-out must be non-decreasing: every month must hold at least as much as the month before. build_out.shift(datetime=1) moves the series forward by one month, and the constraint compares each month against its predecessor in a single vectorized call.

Meeting the energy need#

Multiplying the availability profiles by the installed capacity gives the generation of each resource in each month; broadcasting over datetime and resource handles the whole grid at once. Summing over resource gives the total generation, which must be at least the energy need in every month.

Cost and objective#

Capacity is paid for when it is built, so the cost driver is the capacity added each month — build_out minus the previous month’s value — priced at that month’s unit cost. (In the first month the shift contributes nothing, so the entire initial build-out is charged.) Summing over all months and resources gives the total cost, which becomes the objective to minimize (the default sense of add_objective).

[11]:
build_out = m.add_variables(
    lower=0,
    dims=resource_profiles_xr.dims,
    coords=resource_profiles_xr.coords,
    name="build_out",
)
m.add_constraints(build_out >= build_out.shift(datetime=1), name="build_out_increasing")
/tmp/ipykernel_2254/937858007.py:7: LinopySemanticsWarning: Variable 'build_out' has absent slots (from `mask=` / `.where()` / `.shift()` / `.reindex()`). Under legacy each absent slot contributes 0 to the resulting expression's terms (so `x + y >= 10` reduces to `x >= 10` there). Under v1 the absence propagates through arithmetic instead (`x + y` becomes absent at the slot and the constraint drops).
  Resolve:   wrap with `build_out.fillna(0)` for the legacy behaviour under v1
             (no fix needed if you only use the variable in a constraint LHS alone — `y >= 0` drops the same way in both).
  Opt in:    linopy.options['semantics'] = 'v1'
  Silence:   warnings.filterwarnings('ignore', category=LinopySemanticsWarning)
  m.add_constraints(build_out >= build_out.shift(datetime=1), name="build_out_increasing")
[11]:
Constraint `build_out_increasing` [datetime: 120, resource: 6]:
---------------------------------------------------------------
[2024-01-31 00:00:00, wind_A]: +1 build_out[2024-01-31 00:00:00, wind_A]                                               ≥ -0.0
[2024-01-31 00:00:00, solar_A]: +1 build_out[2024-01-31 00:00:00, solar_A]                                             ≥ -0.0
[2024-01-31 00:00:00, wind_B]: +1 build_out[2024-01-31 00:00:00, wind_B]                                               ≥ -0.0
[2024-01-31 00:00:00, solar_B]: +1 build_out[2024-01-31 00:00:00, solar_B]                                             ≥ -0.0
[2024-01-31 00:00:00, wind_C]: +1 build_out[2024-01-31 00:00:00, wind_C]                                               ≥ -0.0
[2024-01-31 00:00:00, solar_C]: +1 build_out[2024-01-31 00:00:00, solar_C]                                             ≥ -0.0
[2024-02-29 00:00:00, wind_A]: +1 build_out[2024-02-29 00:00:00, wind_A] - 1 build_out[2024-01-31 00:00:00, wind_A]    ≥ -0.0
                ...
[2033-11-30 00:00:00, solar_C]: +1 build_out[2033-11-30 00:00:00, solar_C] - 1 build_out[2033-10-31 00:00:00, solar_C] ≥ -0.0
[2033-12-31 00:00:00, wind_A]: +1 build_out[2033-12-31 00:00:00, wind_A] - 1 build_out[2033-11-30 00:00:00, wind_A]    ≥ -0.0
[2033-12-31 00:00:00, solar_A]: +1 build_out[2033-12-31 00:00:00, solar_A] - 1 build_out[2033-11-30 00:00:00, solar_A] ≥ -0.0
[2033-12-31 00:00:00, wind_B]: +1 build_out[2033-12-31 00:00:00, wind_B] - 1 build_out[2033-11-30 00:00:00, wind_B]    ≥ -0.0
[2033-12-31 00:00:00, solar_B]: +1 build_out[2033-12-31 00:00:00, solar_B] - 1 build_out[2033-11-30 00:00:00, solar_B] ≥ -0.0
[2033-12-31 00:00:00, wind_C]: +1 build_out[2033-12-31 00:00:00, wind_C] - 1 build_out[2033-11-30 00:00:00, wind_C]    ≥ -0.0
[2033-12-31 00:00:00, solar_C]: +1 build_out[2033-12-31 00:00:00, solar_C] - 1 build_out[2033-11-30 00:00:00, solar_C] ≥ -0.0
[12]:
gen = resource_profiles_xr * build_out
total_gen = gen.sum(dim="resource")
m.add_constraints(total_gen >= energy_needs_xr, name="meet_energy_need")
[12]:
Constraint `meet_energy_need` [datetime: 120]:
----------------------------------------------
[2024-01-31 00:00:00]: +0.6112 build_out[2024-01-31 00:00:00, wind_A] - 0.02749 build_out[2024-01-31 00:00:00, solar_A] + 0.2509 build_out[2024-01-31 00:00:00, wind_B] - 0.02102 build_out[2024-01-31 00:00:00, solar_B] + 0.3141 build_out[2024-01-31 00:00:00, wind_C] + 0.03585 build_out[2024-01-31 00:00:00, solar_C] ≥ -0.0
[2024-02-29 00:00:00]: +0.6429 build_out[2024-02-29 00:00:00, wind_A] + 0.109 build_out[2024-02-29 00:00:00, solar_A] + 0.2374 build_out[2024-02-29 00:00:00, wind_B] + 0.2401 build_out[2024-02-29 00:00:00, solar_B] + 0.3239 build_out[2024-02-29 00:00:00, wind_C] + 0.2501 build_out[2024-02-29 00:00:00, solar_C]     ≥ -0.0
[2024-03-31 00:00:00]: +0.4859 build_out[2024-03-31 00:00:00, wind_A] + 0.255 build_out[2024-03-31 00:00:00, solar_A] + 0.2714 build_out[2024-03-31 00:00:00, wind_B] + 0.3011 build_out[2024-03-31 00:00:00, solar_B] + 0.5315 build_out[2024-03-31 00:00:00, wind_C] + 0.4185 build_out[2024-03-31 00:00:00, solar_C]     ≥ -0.0
[2024-04-30 00:00:00]: +0.4231 build_out[2024-04-30 00:00:00, wind_A] + 0.5842 build_out[2024-04-30 00:00:00, solar_A] + 0.4048 build_out[2024-04-30 00:00:00, wind_B] + 0.6578 build_out[2024-04-30 00:00:00, solar_B] + 0.4155 build_out[2024-04-30 00:00:00, wind_C] + 0.5967 build_out[2024-04-30 00:00:00, solar_C]    ≥ -0.0
[2024-05-31 00:00:00]: +0.428 build_out[2024-05-31 00:00:00, wind_A] + 0.7825 build_out[2024-05-31 00:00:00, solar_A] + 0.5631 build_out[2024-05-31 00:00:00, wind_B] + 0.5747 build_out[2024-05-31 00:00:00, solar_B] + 0.3902 build_out[2024-05-31 00:00:00, wind_C] + 0.6521 build_out[2024-05-31 00:00:00, solar_C]     ≥ -0.0
[2024-06-30 00:00:00]: +0.2262 build_out[2024-06-30 00:00:00, wind_A] + 0.995 build_out[2024-06-30 00:00:00, solar_A] + 0.5804 build_out[2024-06-30 00:00:00, wind_B] + 0.9389 build_out[2024-06-30 00:00:00, solar_B] + 0.4771 build_out[2024-06-30 00:00:00, wind_C] + 0.9253 build_out[2024-06-30 00:00:00, solar_C]     ≥ 0.16506350946109716
[2024-07-31 00:00:00]: +0.3046 build_out[2024-07-31 00:00:00, wind_A] + 0.9882 build_out[2024-07-31 00:00:00, solar_A] + 0.6752 build_out[2024-07-31 00:00:00, wind_B] + 0.9205 build_out[2024-07-31 00:00:00, solar_B] + 0.5471 build_out[2024-07-31 00:00:00, wind_C] + 0.8375 build_out[2024-07-31 00:00:00, solar_C]    ≥ 0.6
                ...
[2033-06-30 00:00:00]: +0.2494 build_out[2033-06-30 00:00:00, wind_A] + 0.9619 build_out[2033-06-30 00:00:00, solar_A] + 0.5522 build_out[2033-06-30 00:00:00, wind_B] + 0.9985 build_out[2033-06-30 00:00:00, solar_B] + 0.4841 build_out[2033-06-30 00:00:00, wind_C] + 0.974 build_out[2033-06-30 00:00:00, solar_C]     ≥ 10.9650635094611
[2033-07-31 00:00:00]: +0.3291 build_out[2033-07-31 00:00:00, wind_A] + 1.052 build_out[2033-07-31 00:00:00, solar_A] + 0.6557 build_out[2033-07-31 00:00:00, wind_B] + 0.9252 build_out[2033-07-31 00:00:00, solar_B] + 0.4171 build_out[2033-07-31 00:00:00, wind_C] + 0.8311 build_out[2033-07-31 00:00:00, solar_C]     ≥ 11.400000000000002
[2033-08-31 00:00:00]: +0.2911 build_out[2033-08-31 00:00:00, wind_A] + 0.8516 build_out[2033-08-31 00:00:00, solar_A] + 0.5541 build_out[2033-08-31 00:00:00, wind_B] + 0.9003 build_out[2033-08-31 00:00:00, solar_B] + 0.4281 build_out[2033-08-31 00:00:00, wind_C] + 0.8716 build_out[2033-08-31 00:00:00, solar_C]    ≥ 11.834936490538903
[2033-09-30 00:00:00]: +0.4127 build_out[2033-09-30 00:00:00, wind_A] + 0.5509 build_out[2033-09-30 00:00:00, solar_A] + 0.5896 build_out[2033-09-30 00:00:00, wind_B] + 0.5384 build_out[2033-09-30 00:00:00, solar_B] + 0.3776 build_out[2033-09-30 00:00:00, wind_C] + 0.5582 build_out[2033-09-30 00:00:00, solar_C]    ≥ 13.765063509461093
[2033-10-31 00:00:00]: +0.3214 build_out[2033-10-31 00:00:00, wind_A] + 0.3314 build_out[2033-10-31 00:00:00, solar_A] + 0.4453 build_out[2033-10-31 00:00:00, wind_B] + 0.3801 build_out[2033-10-31 00:00:00, solar_B] + 0.342 build_out[2033-10-31 00:00:00, wind_C] + 0.4634 build_out[2033-10-31 00:00:00, solar_C]     ≥ 16.7
[2033-11-30 00:00:00]: +0.5976 build_out[2033-11-30 00:00:00, wind_A] + 0.2405 build_out[2033-11-30 00:00:00, solar_A] + 0.3778 build_out[2033-11-30 00:00:00, wind_B] + 0.2574 build_out[2033-11-30 00:00:00, solar_B] + 0.3693 build_out[2033-11-30 00:00:00, wind_C] + 0.1193 build_out[2033-11-30 00:00:00, solar_C]    ≥ 18.29519052838329
[2033-12-31 00:00:00]: +0.5195 build_out[2033-12-31 00:00:00, wind_A] + 0.1507 build_out[2033-12-31 00:00:00, solar_A] + 0.3259 build_out[2033-12-31 00:00:00, wind_B] + 0.1729 build_out[2033-12-31 00:00:00, solar_B] + 0.4061 build_out[2033-12-31 00:00:00, wind_C] + 0.0608 build_out[2033-12-31 00:00:00, solar_C]    ≥ 16.5650635094611
[13]:
build_month = build_out - build_out.shift(datetime=1)
cost = build_month * unit_cost_profiles_xr
total_cost = cost.sum()
total_cost
/tmp/ipykernel_2254/918857985.py:1: LinopySemanticsWarning: Variable 'build_out' has absent slots (from `mask=` / `.where()` / `.shift()` / `.reindex()`). Under legacy each absent slot contributes 0 to the resulting expression's terms (so `x + y >= 10` reduces to `x >= 10` there). Under v1 the absence propagates through arithmetic instead (`x + y` becomes absent at the slot and the constraint drops).
  Resolve:   wrap with `build_out.fillna(0)` for the legacy behaviour under v1
             (no fix needed if you only use the variable in a constraint LHS alone — `y >= 0` drops the same way in both).
  Opt in:    linopy.options['semantics'] = 'v1'
  Silence:   warnings.filterwarnings('ignore', category=LinopySemanticsWarning)
  build_month = build_out - build_out.shift(datetime=1)
/home/docs/checkouts/readthedocs.org/user_builds/linopy/envs/latest/lib/python3.12/site-packages/linopy/expressions.py:1243: LinopySemanticsWarning: Coordinate order mismatch in this operator's constant operand: the same labels in a different order were reindexed by label by legacy. Under v1 this raises ValueError (§8).
  Dim:       'resource': left=['wind_A', 'solar_A', 'wind_B', 'solar_B', 'wind_C', 'solar_C'], right=['wind_A', 'wind_B', 'wind_C', 'solar_A', 'solar_B', 'solar_C']
  Resolve:   `.sel(...)` / `.reindex(...)` / `.sortby(...)` to align
             or pass an explicit `join=` argument.
  Opt in:    linopy.options['semantics'] = 'v1'
  Silence:   warnings.filterwarnings('ignore', category=LinopySemanticsWarning)
  self_const, factor, needs_reindex = self._broadcast_and_align(
[13]:
LinearExpression
----------------
+1 build_out[2024-01-31 00:00:00, wind_A] + 0.8 build_out[2024-01-31 00:00:00, solar_A] + 1.2 build_out[2024-01-31 00:00:00, wind_B] ... -0.6067 build_out[2033-11-30 00:00:00, wind_C] + 0.4964 build_out[2033-12-31 00:00:00, solar_C] - 0.4964 build_out[2033-11-30 00:00:00, solar_C]
[14]:
m.add_objective(total_cost)
m
[14]:
Linopy LP model
===============

Variables:
----------
 * build_out (datetime, resource)

Expressions:
------------
<empty>

Constraints:
------------
 * build_out_increasing (datetime, resource)
 * meet_energy_need (datetime)

Status:
-------
initialized

The solution holds the optimal build_out for every resource and month. Reshaping it into a wide table lets us plot how each resource’s installed capacity grows over time.

Solve#

[15]:
m.solve()
Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-4505xesx.lp
Reading time = 0.00 seconds
obj: 840 rows, 720 columns, 2154 nonzeros
Gurobi Optimizer version 13.0.3 build v13.0.3rc0 (linux64 - "Ubuntu 24.04 LTS")

CPU model: AMD EPYC 7R13 Processor, instruction set [SSE2|AVX|AVX2]
Thread count: 1 physical cores, 2 logical processors, using up to 2 threads

Optimize a model with 840 rows, 720 columns and 2154 nonzeros (Min)
Model fingerprint: 0x4ea12f79
Model has 720 linear objective coefficients
Coefficient statistics:
  Matrix range     [2e-03, 1e+00]
  Objective range  [2e-03, 7e-01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e-01, 2e+01]

Presolve removed 689 rows and 112 columns
Presolve time: 0.01s
Presolved: 151 rows, 608 columns, 4278 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0   -1.2611601e+31   2.007125e+33   2.522320e+01      0s
     189    3.4917508e+01   0.000000e+00   0.000000e+00      0s

Solved in 189 iterations and 0.01 seconds (0.00 work units)
Optimal objective  3.491750819e+01
[15]:
('ok', 'optimal')
[16]:
sol = m.solution
[17]:
sol_df = (
    sol.to_dataframe()
    .reset_index()
    .pivot_table(values="build_out", columns="resource", index="datetime")
)
sol_df.plot()
plt.grid()
plt.title("Least Cost Energy Portfolio")
[17]:
Text(0.5, 1.0, 'Least Cost Energy Portfolio')
_images/energy-resource-build_30_1.svg

The plot shows the cumulative installed capacity of each resource over time. The lines only ever rise — capacity is added, never retired — and the optimizer meets each month’s energy need with the cheapest combination of resources.

Conclusion#

By indexing the variable and constraints over datetime and resource, linopy expresses the whole capacity-expansion problem — hundreds of variables and constraints — in a handful of vectorized lines that mirror how the data is laid out in xarray.