Linear Equations

Linear Equations#

Motivation#

Linear algebra is a fundamental tool in scientific computing. Linear systems of equations is one of the most common problems in scientific computing. We will use this module to introduce the basic concept, algorithms and tools for solving linear equations in Python.

Basic Concepts#

If you are given a set of linear equations, you can write it in a matrix form as

\[ A x = b \]

where \(A\) is a matrix, \(x\) is a vector, and \(b\) is a vector. \(A\) is also called the coefficient matrix, and \(b\) is called the right-hand side vector. In linear algebra, \(A\) and \(b\) are often given, and we need to find \(x\).

If the number of equations is equal to the number of unknowns, i.e., the size of \(A\) is \(n\times n\), then \(A\) is called a square matrix. If \(A\) is a square matrix and the determinant of \(A\) is not zero (or say, the rows or columns are linearly independent), then \(A\) is called a nonsingular matrix, and the solution \(x\) is unique.

Matrix Formulation#

Let’s start with a simple example:

\[\begin{split} \begin{bmatrix} 2 & 1 & 4 \\ 1 & 2 & 3 \\ 3 & 1 & 2 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix} \end{split}\]

In Python, the first thing to know is to use numpy to work with matrices and solve the equations.

import numpy as np


A = np.array([[2, 1, 4], [1, 2, 3], [3, 1, 2]])
b = np.array([1, 2, 3])

Some notes:

  • It is very common to do import numpy as np at the beginning of a Python script or a Jupyter notebook.

  • If you are familiar with MATLAB, you may recognize that np.array([]) is the same as [ ] in MATLAB. It can be used to create an array (1D) or an array of arrays (2D).

Solution by Inversion#

The solution of the linear equations \(A x = b\) is, mathematically, \(x = A^{-1} b\). It can be implemented in Python as follows:

Ainv = np.linalg.inv(A)

x = Ainv @ b

print(x)

[ 0.90909091  1.36363636 -0.54545455]

Additional notes are:

  • The np.linalg.inv function is used to compute the inverse of a matrix.

  • The @ operator is used for matrix multiplication in Python.

A quick cheatsheet can be found here for MATLAB users to get started with Python.

But it is never a good idea to use matrix inverse to solve the linear equations, because it can be

  • numerically unstable

  • computationally expensive

For interested readers, here’s a good writeup on why one should never invert a matrix for Ax=b.

Non-Inversion Method#

Instead, we can use the np.linalg.solve function to solve the linear equations.

x = np.linalg.solve(A, b)

print(x)

[ 0.90909091  1.36363636 -0.54545455]

np.linalg.solve uses the LU decomposition to solve the equations. It needs to be understood that the actual decomposition is performed by an underlying library, LAPACK, written in C/Fortran. This is a good example of how Python can be used as a “glue” language to combine different libraries and offer an high-level interface.

LU Decomposition#

LU decomposition (also called LU factorization) is a way to turn a square matrix \(A\) into a lower triangular matrix \(L\) and an upper triangular matrix \(U\), so that \(A = PLU\), where P is a permutation matrix, and L and U are lower and upper triangular matrices.

SciPy allows one to decompose a matrix into LU with scipy.linalg.lu and reapply the decomposition to solve different right-hand side vectors.

from scipy.linalg import lu, lu_factor, lu_solve

lu, piv = lu_factor(A)
# `piv` is the pivot indices; you can read more here:
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_factor.html

x = lu_solve((lu, piv), b)

print(x)

[ 0.90909091  1.36363636 -0.54545455]

If the A matrix stays the same but b changes each time, you can reuse (lu, piv) by calling lu_solve multiple times. This will save the computational cost of repeated LU decomposition, which is the most expensive part in solving linear equations.

Sparse Linear Equations#

A set of equations is called sparse if the \(A\) matrix has many zero elements. When using sparse matrix types to store \(A\), the linear equations can be solved with better efficiency and less memory usage.

Sparse matrices are common in power system applications. As is typical, a substation represented as a bus is only connected to a few other buses. If we write out the nodal equation \(I = Y V\), the admittance matrix \(Y\) is sparse.

Let’s use the numbers from the previous example but create a sparse matrix with scipy.

A = np.array([[2, 1, 4], [1, 2, 3], [3, 1, 2]])
b = np.array([1, 2, 3])

from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve

A_sparse = csr_matrix(A)

x = spsolve(A_sparse, b)

print(x)

[ 0.90909091  1.36363636 -0.54545455]

A few notes are:

  • The sparse matrix is created by csr_matrix from a dense matrix \(A\). This is not the most efficient way if the sparse matrix can be created without creating a full matrix. Here we use the dense matrix to illustrate the solution.

  • scipy.sparse.linalg.spsolve is used to solve the sparse linear equations. The solution is the same as the dense matrix.

If we inspect A_sparse, we can see that it is a sparse matrix.

print(A_sparse)

<Compressed Sparse Row sparse matrix of dtype 'int64'
	with 9 stored elements and shape (3, 3)>
  Coords	Values
  (0, 0)	2
  (0, 1)	1
  (0, 2)	4
  (1, 0)	1
  (1, 1)	2
  (1, 2)	3
  (2, 0)	3
  (2, 1)	1
  (2, 2)	2

SciPy supports many types of sparse matrices, including csr_matrix, csc_matrix, coo_matrix, dia_matrix, dok_matrix, lil_matrix. A summary can be found here.

The differences are in the underlying storage for representing the matrix on a computer. For example, coo_matrix stores the matrix in a coordinate format, which is efficient for creating a sparse matrix from a dense matrix. Similarly, dok_matrix stores the matrix in a dictionary, which is also efficient for increamentally building a matrix. These two formats are not suitable for numerical operations because of the difficulty in accessing the elements in any desired order. Specifically, for example, to access (1, 1) in a coo_matrix, we need to search through all the stored elements.

Working with Sparse Matrix#

As discussed, creating a sparse matrix from dense is not efficient, because when creating the dense matrix, heavy memory allocations are performed. Instead, we need to understand the original problem well and create the sparse matrix directly.

Consider a four-bus power system with the line parameters given in the table below. R and X are the series resistance and reactance. B is the shunt susceptance which should be divided by 2 for each terminal.

From

To

R

X

B

0

1

0.01

0.05

0.02

0

2

0.02

0.06

0.03

1

2

0.03

0.07

0.04

2

3

0.04

0.08

0.05

For the sake of simplicity, we start the bus index with 0, which is consistent with the Python index.

We build the admittance matrix using a dok_matrix format.

from scipy.sparse import dok_matrix
import numpy as np

r = np.array([0.01, 0.02, 0.03, 0.04])
x = np.array([0.05, 0.06, 0.07, 0.08])
b = np.array([0.02, 0.03, 0.04, 0.05])

fr = np.array([0, 0, 1, 2])
to = np.array([1, 2, 2, 3])

nlines = 4

# create an empty matrix with 4 rows and 4 columns
Ybus = dok_matrix((4, 4), dtype=complex)

Internally, the dok_matrix is a dictionary that maps the row index and column index to the value. It is useful to build the matrix incrementally, because you will be doing random access to find the elements to update. Based on hash table, dictionary access is on average O(1).

# follow the definition of admittance and insert the elements

for i in range(nlines):
    Ybus[fr[i], fr[i]] += 1 / (r[i] + 1j * x[i]) + b[i] / 2
    Ybus[to[i], to[i]] += 1 / (r[i] + 1j * x[i]) + b[i] / 2
    Ybus[fr[i], to[i]] -= 1 / (r[i] + 1j * x[i])
    Ybus[to[i], fr[i]] -= 1 / (r[i] + 1j * x[i])

print(Ybus)

<Dictionary Of Keys sparse matrix of dtype 'complex128'
	with 12 stored elements and shape (4, 4)>
  Coords	Values
  (0, 0)	(8.871153846153845-34.230769230769226j)
  (1, 1)	(9.048567639257293-31.299734748010607j)
  (0, 1)	(-3.846153846153846+19.23076923076923j)
  (1, 0)	(-3.846153846153846+19.23076923076923j)
  (2, 2)	(15.232413793103449-37.06896551724138j)
  (0, 2)	(-5.000000000000001+15j)
  (2, 0)	(-5.000000000000001+15j)
  (1, 2)	(-5.1724137931034475+12.068965517241379j)
  (2, 1)	(-5.1724137931034475+12.068965517241379j)
  (3, 3)	(5.025-10j)
  (2, 3)	(-5+10j)
  (3, 2)	(-5+10j)

One will then convert the dok_matrix to a coo_matrix to perform the matrix-vector multiplication.

Ybus_coo = Ybus.tocsr()

print(Ybus_coo)

<Compressed Sparse Row sparse matrix of dtype 'complex128'
	with 12 stored elements and shape (4, 4)>
  Coords	Values
  (0, 0)	(8.871153846153845-34.230769230769226j)
  (0, 1)	(-3.846153846153846+19.23076923076923j)
  (0, 2)	(-5.000000000000001+15j)
  (1, 1)	(9.048567639257293-31.299734748010607j)
  (1, 0)	(-3.846153846153846+19.23076923076923j)
  (1, 2)	(-5.1724137931034475+12.068965517241379j)
  (2, 2)	(15.232413793103449-37.06896551724138j)
  (2, 0)	(-5.000000000000001+15j)
  (2, 1)	(-5.1724137931034475+12.068965517241379j)
  (2, 3)	(-5+10j)
  (3, 3)	(5.025-10j)
  (3, 2)	(-5+10j)

Each sparse matrix format in SciPy has supports different ways to construct the matrix. Please refer to the documentation for the specifics.

Question: What other ways can you think of to build the admittance matrix? How about not using the for loop?

Alternative Sparse Solvers#

The default underlying solver for scipy.sparse.linalg.spsolve is SuperLU. You can install the Python package scikit-sparse to use a different solver, UMFPACK, which is the default solver for MATLAB. The two solvers are inherently different. One needs to benchmark the solvers for the specific problem to find out the more efficient one.

Below is an example of benchmarking the two solvers.

# make a large random sparse matrix
from scipy.sparse import rand
import numpy as np
A = rand(10000, 10000, density=0.0002, format='dok')

Question: Test how long it takes if the format is changed to csr.

Testing on a power system case#

We first download the case file from MATPOWER.

import urllib.request
import os
import tempfile

# URL of the file
url = "https://raw.githubusercontent.com/MATPOWER/matpower/refs/heads/master/data/case9241pegase.m"

# Create a temporary file
with tempfile.NamedTemporaryFile(delete=False, suffix='.m') as temp_file:
    # Download the file
    urllib.request.urlretrieve(url, temp_file.name)

    print(f"File downloaded and saved to {temp_file.name}")

    # Verify if the file exists
    if os.path.exists(temp_file.name):
        print(f"File size: {os.path.getsize(temp_file.name)} bytes")
    else:
        print("Download failed")

    # The file and temporary directory will be automatically deleted
    # when we exit the 'with' block

File downloaded and saved to /tmp/tmpgjvea826.m
File size: 1530051 bytes

Next, load the case into ANDES and obtain the Jacobian matrix for power flow at the first iteration.

import andes
from andes.linsolvers.scipy import spmatrix_to_csc


ss = andes.load(temp_file.name, no_output=True)
ss.PFlow.init()
ss.PFlow.nr_step()

# get the Jacobian matrix
J = spmatrix_to_csc(ss.dae.gy)

# get the right-hand side vector
b = np.array(ss.dae.g)

We can verify the size of the Jacobian matrix:

J

<Compressed Sparse Column sparse matrix of dtype 'float64'
	with 154958 stored elements and shape (19928, 19928)>

# Solve the linear equations with SuperLU

from scipy.sparse.linalg import spsolve
%timeit -n 20 -r 5 spsolve(J, b, use_umfpack=False)

27 ms ± 560 μs per loop (mean ± std. dev. of 5 runs, 20 loops each)

In the above, we passed two arguments to %timeit to run the code in a loop of 20 times and repeat 5 loops. This is to get a more accurate measurement of the execution time. This manual adjustment is useful when the default %timeit takes too long to run.

Next, let’s solve it with UMFPACK. Make sure you have installed scikit-umfpack before running the following code.

from scikits.umfpack import spsolve as umfpack_spsolve


%timeit -n 20 -r 5 umfpack_spsolve(J, b)

36.6 ms ± 1.13 ms per loop (mean ± std. dev. of 5 runs, 20 loops each)

The result shows that UMFPACK takes more time to run.

This is not to say that UMFPACK is always slower. The two solvers are designed for sparse matrices with different underlying structures. Again, the performance needs to be benchmarked for specific cases.

Alternative Sparse Libraries#

In the end, we would like to mention that NumPy and SciPy are not the only libraries for scientific computing in Python. Without doubt, they are the most popular ones, which means they have been tested over a wide range of cases and are well-documented. Users are less likely to run into bugs or weird corner cases, even though they are supported by the open-source community. In addition, the ecosystem around NumPy and SciPy are very rich, meaning that there are other packages built on top of them or interoperate with them v.

There are other libraries for matrix computation, such as CVXOPT and its fork KVXOPT. As its name suggests, CVXOPT is a package designed for convex optimization. The authors of CVXOPT also coded up dense and sparse matrix types and related linear algebra operations.

Like SciPy, CVXOPT does not develop its own sparse linear solvers but interfaces with the existing, highly optimized ones. Over the years, many efficient solvers have been developed, such as SuperLU, UMFPACK, KLU, among many others. CVXOPT interfaces its sparse matrix type with UMFPACK, and KVXOPT further extends it to support KLU, which is “Kent Clark LU”. KLU is worth mentioning because it shows superior performance for electrical circuit problems.

In ANDES, the Jacobian matrix is of the type KVXOPT.spmatrix. It can work with both the UMFPACK and KLU solvers.

Let’s prepare the data, J and b matrices for the KLU solver. Note that the matrix b needs to be converted to the matrix type to be supported by KVXOPT. Fortunately, under the hood, matrix is interoperable with NumPy arrays.

from kvxopt import klu, umfpack, matrix


J = ss.dae.gy

# convert b to KVXOPT.matrix type
b = matrix(ss.dae.g)

Let’s first test it with UMFPACK:

%%timeit
b_new = matrix(b)
umfpack.linsolve(J, b_new)

37.4 ms ± 578 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In the above, %%timeit is a magic command to time the cell, which can be multiple lines. We have to create a copy of b every time because, by the design of the solver, the output is written to the b matrix. Creating a copy of b ensures that we are solving the same problem each time.

If you compare the UMFPACK time via SciPy and KVXOPT, you will find that they are very close. This is largely because both of the wrappers in Python are efficient, and most of the time is spend on the underlying solver, which is the same via SciPy and KVXOPT.

Next, let’s test it with KLU:

%%timeit
b_new = matrix(b)
klu.linsolve(J, b_new)

10 ms ± 57.9 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

KLU only takes about 1/3 of the time of UMFPACK. Again, this is not to say that KLU is always better. They suit different types of problems, and the way for the user to choose the right solver is to be practical: benchmark them for the specific problems.

Note that if you are creating the sparse matrix with SciPy, there is an overhead to convert it to the KVXOPT.spmatrix type.

Exercise: Please read the documentation for CVXOPT/KVXOPT and write a function to convert a SciPy sparse matrix to a KVXOPT.spmatrix. Which SciPy sparse matrix will you choose?

In addition, sparse LU factorization results can be reused for different right-hand sides, just like in the dense canse. This will be discussed in the section for acceleration techniques.

Further Reading#

There are iterative approaches for solving linear equations. Iterative solvers have not seen widespread applications in power system tools, so we leave the reading to interested readers.

Summary#

  • Fundamentally, we can solve linear equations \(A x = b\) in Python with numpy.linalg.solve or scipy.sparse.linalg.spsolve.

  • For large-scale sparse problems, the sparse matrix solvers are more efficient.

  • We can benchmark the solvers for specific problems to find the most efficient one.

  • The choice of sparse matrix libraries can also impact the performance.

Contents