Integrated Energy System Planning#

The application domain of Integrated Energy System Planning (IESP) is a field that describes the energy system in mathematical models typically formulated as optimization problems. Solving energy planning tasks is essential in the transition to climate neutrality of regional and global energy systems, providing decision support to policymakers by gaining insights into complex interactions and dynamics of increasingly integrated energy systems, e.g., see Schmitz et al., 2023, Böttger and Härtel 2022, Craig et al. 2022, Frischmuth et al. 2022.

The planning problems corresponding to large-scale energy systems (easily up to 100 million decision variables) are typically formulated as continuous linear programming (LP) problems. Note that there are many non-linear and non-convex properties in energy systems. However, especially in long-term analyses with a large geographical scope, modeling applications apply simplifications, abstractions, and approximations. They resort to (convex) linear programming formulations to reduce the computational burden.

[2]:

# Uncomment the following line to install dprox if your are in the online google colab notebook
# !pip install dprox

A reduced IESP problem instance#

We consider an exemplary LP problem instance representing an IESP problem. It models a (reduced) multi-period capacity expansion and system operation planning problem instance which helps to develop transformation pathways of integrated energy systems.

Importing the required packages#

[2]:

import os
import sys
from scipy import io
import scipy.sparse
import numpy as np
import dprox.utils.huggingface as hf

Load problem instance#

Let us load a small IESP problem example in standard LP format. Note that the constraint matrix is stored in a sparse format and distinguishes equality (eq) and inequality (ub) constraints.

[4]:

def load_lp_componentes_from_mat(filename):
    """
    Loads the linear program components from a .mat file.
    """
    model_components = io.loadmat(filename)
    A_ub = model_components["A_ub"]
    b_ub = model_components["b_ub"][0, :]
    A_eq = model_components["A_eq"]
    b_eq = model_components["b_eq"][0, :]
    c = model_components["c"][0, :]
    x_lb = model_components["x_lb"][0, :]
    x_ub = model_components["x_ub"][0, :]

    return (A_ub, b_ub, A_eq, b_eq, c, x_lb, x_ub)

# Load LP components
A_ub, b_ub, A_eq, b_eq, c, x_lb, x_ub = load_lp_componentes_from_mat(hf.load_path("examples/iesp/esm_instance.mat"))

Problem size and coefficient statistics#

Next, we report the problem size and coefficient statistics for the energy system planning problem. The problem size is reported in terms of the number of variables and constraints. The coefficient statistics are reported by the minimum and maximum values of non-zero coefficients in the corresponding linear problem matrices and vectors.

[17]:

# Report problem size
print(f"Number of decision variables: {A_eq.shape[1]}")
print(f"Number of linear equality constraints: {A_eq.shape[0]}")
print(f"Number of linear inequality constraints (w/o bound constraints): {A_ub.shape[0]}\n")

# Report coefficient statistics
print("Coefficient statistics (positive non-zero) of model components:")
print(f"c \t: [{np.min(np.abs(c[c != 0])):3.2e}, {np.max(np.abs(c[c != 0])):3.2e}]")
print(f"A_eq \t: [{np.min(np.abs(A_eq[A_eq != 0])):3.2e}, {np.max(np.abs(A_eq[A_eq != 0])):3.2e}]")
print(f"A_ub \t: [{np.min(np.abs(A_ub[A_ub != 0])):3.2e}, {np.max(np.abs(A_ub[A_ub != 0])):3.2e}]")
print(f"b_eq \t: [{np.min(np.abs(b_eq[b_eq != 0])):3.2e}, {np.max(np.abs(b_eq[b_eq != 0])):3.2e}]")
print(f"b_ub \t: [{np.min(np.abs(b_ub[b_ub != 0])):3.2e}, {np.max(np.abs(b_ub[b_ub != 0])):3.2e}]")
print(f"x_ub \t: [{np.min(np.abs(x_ub[x_ub != 0])):3.2e}, {np.max(np.abs(x_ub[x_ub != 0])):3.2e}]")

Number of decision variables: 324150
Number of linear equality constraints: 52558
Number of linear inequality constraints (w/o bound constraints): 306604

Coefficient statistics (positive non-zero) of model components:
c       : [1.00e-08, 4.38e-01]
A_eq    : [9.59e-01, 1.07e+00]
A_ub    : [1.01e-06, 6.00e+00]
b_eq    : [1.00e+00, 1.96e+05]
b_ub    : [4.80e+03, 2.00e+08]
x_ub    : [inf, inf]

Visualize the sparsity pattern of the linear problem matrices#

Next, we visualize the sparsity pattern of the linear problem matrices. This is useful to get an idea of the structure of the problem and to check if the problem is formulated correctly.

[ ]:

import matplotlib.pyplot as plt
def plot_sparsity_pattern(A, name):
    plt.spy(A, markersize=0.25, aspect="auto")
    plt.title(f"Sparsity pattern of the {name} matrix (sparsity={A.getnnz() / (A.shape[0] * A.shape[1]):.2e})")

plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plot_sparsity_pattern(A_ub, "A_ub")
plt.subplot(1, 2, 2)
plot_sparsity_pattern(A_eq, "A_eq")

Note that the energy system planning instance poses an ill-conditioned and sparse linear programming problem. Due to the different commodities, technologies, and actors involved, this is a typical characteristic of energy system planning problem instances.

Solving the IESP problem with \(\nabla\)-Prox#

Preparation#

First, we have to import the necessary libraries, configure the CUDA device and set necessary algorithm parameters.

[7]:

os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"
import time
import random
import torch
import dprox as dp

# Set device
if torch.cuda.device_count() > 0:
    device = torch.device('cuda')
else:
    device = torch.device('cpu')

# Set seed for reproducibility
random.seed(2023)
torch.manual_seed(2023)
np.random.seed(2023)

# Set algorithm parameters
max_iters = 100000
abs_tol = 5e-2
rel_tol = 1e-10
rho = 1e-2
dtype = torch.float64
norm_ord = float('inf')

Solve the problem with \(\nabla\)-Prox#

[ ]:

# Initialize problem and solver
start = time.time()
lp_problem = dp.algo.lp.LPProblem(c, A_ub, b_ub, A_eq, b_eq, x_lb, x_ub, norm_ord=norm_ord, dtype=dtype, sparse=True, device=device)
lp_admm = dp.algo.lp.LPSolverADMM(rho=rho, problem_scale=None, abstol=abs_tol, reltol=rel_tol, max_iters=max_iters, dtype=dtype).to(device)

# Solve problem
with torch.no_grad():
    lp_admm.eval()
    x, history, res = lp_admm.solve(lp_problem, residual_balance=True)

x = x.cpu().numpy().squeeze()

print(res[0])

print(f'x >= 0 err: {x.min():.2e}')
print(f'Ax = b err: {np.abs(b_eq - A_eq @ x).max():.2e}')
print(f'Ax <= b err: {np.min(b_ub - A_ub @ x):.2e}')

torch.cuda.current_stream().synchronize()
end = time.time()
print(f"Time elapsed (LPADMM): {end-start:.4f}s\n")

Solving with dprox on a NVIDIA A100 GPU results in the following output:

tensor([[83432.3454]], device='cuda:0', dtype=torch.float64)
x >= 0 err: -1.90e-02
Ax = b err: 1.33e-02
Ax <= b err: -5.98e-02
Time elapsed (LPADMM): 247.9268s

Comparison with established CPU solvers#

To compare the results and convergence behavior with established solvers, we use the proprietary Gurobi Optimizer and open-source HiGHS solvers.

Gurobi#

First, we solve the problem instance with the Gurobi solver via the gurobipy interface. Please note that this step requires a valid Gurobi license as it is a proprietary solver.

[ ]:

import gurobipy as gp
# Requires a license

m = gp.Model()
x = m.addMVar(c.shape[0], lb=x_lb, ub=x_ub)
m.setObjective(c @ x)
m.addConstr(A_eq @ x == b_eq)
m.addConstr(A_ub @ x <= b_ub)

m.params.Method = 2 # Interior Point Method (Barrier + crossover)

m.optimize()

Solving with Gurobi Optimizer version 10.0.2 build v10.0.2rc0 (win64)

CPU model: 12th Gen Intel(R) Core(TM) i7-1280P
Thread count: 14 physical cores, 20 logical processors, using up to 20 threads

[...]

Solved in 292328 iterations and 85.80 seconds (355.13 work units)
Optimal objective  8.341841135e+04

HiGHS#

We use the highspy interface to solve the LP problem with HiGHS. The highspy interface is a python wrapper for the HiGHS solver.

[ ]:

# !pip install highspy

[ ]:

import highspy
import scipy.sparse

h = highspy.Highs()
h.addVars(len(c), x_lb, x_ub)
h.changeColsCost(len(c), np.arange(len(c), dtype=np.int32), c)

A = scipy.sparse.vstack((A_eq, A_ub))
b = np.hstack((b_eq, b_ub))

num_cons = A.shape[0]
A = A.tocsr()

lower=np.hstack((b_eq, -np.inf * np.ones_like(b_ub)))
upper=np.hstack((b_eq, b_ub))

h.addRows(num_cons, lower, upper, A.nnz, A.indptr, A.indices, A.data)

lp = h.getLp()
h.passModel(lp)
h.setOptionValue("log_to_console", True)
h.setOptionValue("solver", "ipm")

start = time.time()

h.run()

end = time.time()
print(f"Time elapsed: {end - start:.2f}s")

info = h.getInfo()
model_status = h.getModelStatus()
print('Model status = ', h.modelStatusToString(model_status))
print('Optimal objective = ', info.objective_function_value)
print('Iteration count = ', info.simplex_iteration_count)
print('Primal solution status = ', h.solutionStatusToString(info.primal_solution_status))
print('Dual solution status = ', h.solutionStatusToString(info.dual_solution_status))
print('Basis validity = ', h.basisValidityToString(info.basis_validity))

Time elapsed: 3555.89s
Model status =  Optimal
Optimal objective =  83418.40550119097
Iteration count =  57338
Primal solution status =  Feasible
Dual solution status =  Feasible
Basis validity =  Valid

Summary#

For the example IESP problem, we see that \(\nabla\)-Prox can solve the problem with solutions times comparable to Gurobi and faster than HiGHS with accuracy close to the results of Gurobi and HiGHS.

Since \(\nabla\)-Prox only relies on first-order methods, it offers substantial memory savings, which is key for large-scale problems as memory use presents a challenge for established CPU-based solvers and their second-order methods.