Linear Algebra in Engineering: Equilibrium in Linear Systems (Part 1 of 7)

Preface
This series is aimed at providing tools for an electrical engineer to analyze data and solve problems in design. The focus is on applying linear algebra to systems of equations or large sets of matrix data.

Introduction
This article will demonstrate the use of matrix algebra to solve for equilibrium in systems. This is common in network flow, economics and electrical circuits (current and voltage analysis). We will apply exact and least squares solutions.

If you are not familiar with linear algebra or need a brush up I recommend Linear Algebra and its Applications by Lay, Lay and McDonald. It provides an excellent review of theory and applications.

This also assumes you are familiar with Python or can stumble your way through it.

The data and code are available here: linear1.py and linear1.xlsx.

Concepts
Linear Equation: an equation that can be written in the form b = a₁x₁ + a₂x₂ + ... + a_nx_n.

Linear System: A set of linear equations involving the same variables.

Solution: a set of values for the set of variables in a linear system that make the system a true statement.

Linear systems are either consistent (having a solution) or inconsistent (having no solution). Linear systems with a solution are either unique (1 solution) or have a free variable indicating infinitely many solutions.

Matrix Equation: The above definition of a linear equation uses a linear combination to create a vector equation. The same equation can be transliterated to matrix notation (used by computers) in the form Ax = b.

Invertible Matrix Theorem: from Lay, Lay, McDonald.
Let A be a square nxn matrix. Then the following statements are equivalent. That is, for a given A, the statements are either all true or all false.

1. A is an invertible matrix.
2. A is row equivalent to the nxn identity matrix.
3. A has n pivot positions.
4. The equation Ax=0 has only the trivial solution.
5. The columns of A form a linearly independent set.
6. The linear transformation x->Ax is one-to-one.
7. The equation Ax=b has at least one solution for each b in Rⁿ.
8. The columns of A span Rⁿ.
9. The linear transformation x->Ax maps Rⁿ onto Rⁿ.
10. There is an nxn matrix C such that CA=I.
11. There is an nxn matrix D such that AD=I.
12. A^T is an invertible matrix.
13. The columns of A form a basis of Rⁿ.
14. Col A=Rⁿ.
15. dim Col A=n
16. rank A=n
17. Nul A={0}
18. dim Nul A=0
19. The number 0 is not an eigenvalue of A.
20. The determinant of A is not zero.

Importing Your Data Set
I will use the Anaconda package suite with Python, numpy and the Spyder IDE for mathematical analysis. It is cross platform, free and open source. It is also easy to import data from Excel once you have the code snippet.

I usually have the following boilerplate code in my scripts:

import sys
import xlrd
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d


np.set_printoptions(threshold=np.nan)
print sys.version
print __name__


The system in this example is traffic flow for 4 intersections in a city grid. All ingress/egres traffic is known and is pictured below.

The data above has been arranged into matrix notation and stored in the excel file with the coefficients in the first five columns and the dependent variables in the sixth (last) column.


filename='linear_1.xlsx'
print
print 'Opening',filename
ss=xlrd.open_workbook(filename,on_demand=True)
sheet_index=0
print ' Worksheet',ss.sheet_names()[sheet_index]
sh=ss.sheet_by_index(sheet_index)


print ' Reading data values'
spec=[]
for row_index in range(0,sh.nrows):
spec.append([sh.cell_value(row_index,0),sh.cell_value(row_index,1),sh.cell_value(row_index,2),sh.cell_value(row_index,3),sh.cell_value(row_index,4)])
spec=np.array(spec)
print ' Design spec array size:',len(spec)


b=[]
for row_index in range(0,sh.nrows):
b.append(sh.cell_value(row_index,5))
b=np.array(b)
print ' Design b array size: ',len(b)


Set filename to your Excel filename and set sheet_index to the tab you're reading from (counting starts from zero). Set range(x,y) and {a,b,c,...,n} to the array you are reading in from the sheet. range(0,sh.nrows) is a handy expression to use here.

The last bit of code converts the array to a numpy array which will be necessary for the code below.

The whole operation has fairly good performance.

Solving the System for Equilibrium
Now that the test results and input variables have been imported let's solve the system. We want to find the equilibrium price.

I have found that np.linalg.solve() is not always the best solution because it cannot handle systems that are underdetermined, overdetermined or linearly dependent. A better alternative seems to be the np.linalg.lstsq() function.

Both functions take the form function(a,b) where a and b are from ax=b.


#eq=np.linalg.solve(spec,b)
eq=np.linalg.lstsq(spec,b)
print
print eq


This is the output:

Opening linear_1.xlsx
Worksheet Sheet1
Reading data values
Design spec array size: 6
Design b array size: 6


(array([ 375., 425., 400., 275., 225.]), array([], dtype=float64), 4, array([ 2.11688316e+00, 1.63978419e+00, 1.41421356e+00,
9.10995891e-01, 8.65181300e-17]))


The first array is the compute result. DO NOT take that result for granted. You must examine the second array for anything unusual. The most relevant thing to note is that the fifth value is close to zero. This is a strong indication of a free (dependent) variable.

Adding a sixth row to the data set that hard sets the fifth variable to 9 results in this output:

Opening linear_1.xlsx
Worksheet Sheet1
Reading data values
Design spec array size: 6
Design b array size: 6


(array([ 591., 209., 400., 491., 9.]), array([ 4.32956906e-27]), 5, array([ 2.16121248, 1.7660084 , 1.47046231, 0.92978281, 0.42850842]))


An entirely different and this time correct result. We can see the fifth variable hard set to 9. If you plug the x array into each row of coefficients you'll find 100% agreement.

Systems of Systems
It is often desirable to solve more than one system. This can be done by iterating through each array, one at a time. This is okay but I have found that memory requirements can be utterly huge for this method. The method below seems bulky but is exponentially faster and uses less memory (I don't know how much because the last time I tried the iterative method I kept running out of memory after 4-8 hours).

The following pseudo code works well when the set of coefficients is the same and the dependent values are different. It works just as well for different sets of coefficients you just have to build your array appropriately.

I also include a determinant check to verify the array is invertible and therefore solvable.


a=np.array([[some],[set],[of],[coefficients]])
b=np.array([[some],[set],[of],[dependents]])
n=number of coefficients
if np.linalg.det(a) != 0:
csi=np.linalg.solve([a]*n,b)


Next Up
In the next article we will examine trend analysis using polynomial least squares.