Portfolio optimization

This page contains advanced examples for the tno.quantum.problems.portfolio_optimization package. Examples of basic usage can be found in the module documentation.

The package provides code that converts the multi-objective portfolio optimization problem into a QUBO problem. The transformed problem can then be solved using quantum annealing techniques.

Requirements

Install the following dependencies to run the examples below:

pip install tno.quantum.problems.portfolio_optimization~=2.0
pip install matplotlib

Examples

Example 1: Demonstrates how to use different QUBO solvers for the portfolio optimization problem.

"""Demonstrate usage of different QUBO solvers for portfolio optimization."""

import matplotlib.pyplot as plt
import numpy as np

from tno.quantum.optimization.qubo import SolverConfig
from tno.quantum.problems.portfolio_optimization import (
    PortfolioOptimizer,
    plot_points,
)

# Set up penalty coefficients for the constraints
lambdas1 = np.logspace(-16, 1, 25, endpoint=False, base=10.0)
lambdas2 = np.logspace(-16, 1, 25, endpoint=False, base=10.0)
lambdas3 = np.array([1])

# Create portfolio optimization problem
portfolio_optimizer = PortfolioOptimizer("benchmark_dataset")
portfolio_optimizer.add_minimize_hhi(weights=lambdas1)
portfolio_optimizer.add_maximize_roc(formulation=1, weights_roc=lambdas2)
portfolio_optimizer.add_emission_constraint(
    weights=lambdas3,
    emission_now="emis_intens_now",
    emission_future="emis_intens_future",
    name="emission",
)

# Define different solver configurations
solver_configurations: dict[str, SolverConfig] = {
    "Simulated Annealing": SolverConfig(
        name="simulated_annealing_solver",
        options={"num_reads": 100, "num_sweeps": 200},
    ),
    "Random Sampler": SolverConfig(
        name="random_sampler_solver",
        options={"num_reads": 100},
    ),
    "Quantum Annealing": SolverConfig(
        name="d_wave_clique_sampler_solver",
        options={"num_reads": 100, "annealing_time": 100},
    ),
}

# Compute the portfolio optimization experiments
experiment_data = []
for solver_config in solver_configurations.values():
    results = portfolio_optimizer.run(solver_config, verbose=False)
    results_df = results.results_df
    x, y, z = (
        results_df["diff diversification"],
        results_df["diff ROC"],
        results_df["diff emission"],
    )
    experiment_data.append((x, y, z))

# Plot the experiment data
fig, axs = plt.subplots(ncols=3, figsize=(12, 5))
(ax1, ax2, ax3) = axs
for i, ((x, y, z), solver_label) in enumerate(
    zip(experiment_data, solver_configurations)
):
    ax = axs[i]
    scatter_plot = plot_points(x, y, c=z, cmap="viridis", ax=ax, vmin=-35, vmax=-25)
    ax.set_title(solver_label)
    ax.set_xlim(-15, 5)
    ax.set_ylim(-3, 10)
cbar = fig.colorbar(scatter_plot, ax=axs, orientation="horizontal", label="Emission")
plt.show()

Additionally we can make a boxplot of emission values.

# Create boxplot of emission values
data = [z for (_, _, z) in experiment_data]
labels = solver_configurations.keys()
plt.boxplot(data, tick_labels=labels)
plt.ylabel("Emission Values")
plt.title("Emission Box plot")
plt.tight_layout()
plt.show()

Example 2: Demonstrates how to use emission constraint in portfolio optimization.

"""Demonstrate usage of emission constraints."""

import matplotlib.pyplot as plt
import numpy as np

from tno.quantum.optimization.qubo import SolverConfig
from tno.quantum.problems.portfolio_optimization import (
    PortfolioOptimizer,
    plot_points,
)

# Define solver configurations
solver_config = SolverConfig(
    name="simulated_annealing_solver",
    options={"num_reads": 20, "num_sweeps": 200},
)

# Set up penalty coefficients for the constraints
lambdas1 = np.logspace(-16, 1, 25, endpoint=False, base=10.0)
lambdas2 = np.logspace(-16, 1, 25, endpoint=False, base=10.0)
lambdas3 = np.array([1])

# Define basic portfolio optimization problem instance
portfolio_optimizer = PortfolioOptimizer(portfolio_data="benchmark_dataset", k=2)
portfolio_optimizer.add_minimize_hhi(weights=lambdas1)
portfolio_optimizer.add_maximize_roc(formulation=1, weights_roc=lambdas2)

# Solve problem without emission constraint
results_no_constraint = portfolio_optimizer.run(solver_config)
results_no_constraint_df = results_no_constraint.results_df

# Compute difference in emission
portfolio_data = results_no_constraint.portfolio_data
outstanding_now = portfolio_data.get_outstanding_now()
total_outstanding_now = np.sum(outstanding_now)

emission_intensity_now = portfolio_data.get_column("emis_intens_now")
emission_intensity_future = portfolio_data.get_column("emis_intens_future")


def calculate_diff_emission(outstanding_future):
    """Calculate manually "diff emission" values"""
    total_outstanding_future = np.sum(outstanding_future)
    total_relative_emission_now = (
        np.sum(outstanding_now * emission_intensity_now) / total_outstanding_now
    )
    total_relative_emission_future = (
        np.sum(outstanding_future * emission_intensity_future)
        / total_outstanding_future
    )
    return 100 * (total_relative_emission_future / total_relative_emission_now - 1)


results_no_constraint_df["diff emission"] = results_no_constraint_df[
    "outstanding amount"
].apply(calculate_diff_emission)

x1, y1, z1 = (
    results_no_constraint_df["diff diversification"],
    results_no_constraint_df["diff ROC"],
    results_no_constraint_df["diff emission"],
)

# Solve problem with emission constraint
portfolio_optimizer.add_emission_constraint(
    weights=lambdas3,
    emission_now="emis_intens_now",
    emission_future="emis_intens_future",
    name="emission",
)
results_constraint = portfolio_optimizer.run(solver_config)
results_constraint_df = results_constraint.results_df

x2, y2, z2 = (
    results_constraint_df["diff diversification"],
    results_constraint_df["diff ROC"],
    results_constraint_df["diff emission"],
)

# Plot the experiment data
fig, axs = plt.subplots(ncols=2, figsize=(12, 5))
(ax1, ax2) = axs

scatter_plot1 = plot_points(x1, y1, c=z1, cmap="viridis", ax=ax1, vmin=-35, vmax=-25)
scatter_plot2 = plot_points(x2, y2, c=z2, cmap="viridis", ax=ax2, vmin=-35, vmax=-25)
cbar = fig.colorbar(scatter_plot1, ax=axs, orientation="horizontal", label="Emission")
ax1.set_title("Unconstrained problem")
ax2.set_title("Constrained problem")

for ax in axs:
    ax.set_xlim(-15, 5)
    ax.set_ylim(-2, 10)

plt.show()