ChEn-3170: Computational Methods in Chemical Engineering Spring 2020 UMass Lowell; Prof. V. F. de Almeida 23Mar20
$ \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}} \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}} \newcommand{\Mmtrx}{\boldsymbol{\mathsf{M}}} \newcommand{\Imtrx}{\boldsymbol{\mathsf{I}}} \newcommand{\Pmtrx}{\boldsymbol{\mathsf{P}}} \newcommand{\Lmtrx}{\boldsymbol{\mathsf{L}}} \newcommand{\Umtrx}{\boldsymbol{\mathsf{U}}} \newcommand{\Smtrx}{\boldsymbol{\mathsf{S}}} \newcommand{\xvec}{\boldsymbol{\mathsf{x}}} \newcommand{\avec}{\boldsymbol{\mathsf{a}}} \newcommand{\bvec}{\boldsymbol{\mathsf{b}}} \newcommand{\cvec}{\boldsymbol{\mathsf{c}}} \newcommand{\rvec}{\boldsymbol{\mathsf{r}}} \newcommand{\mvec}{\boldsymbol{\mathsf{m}}} \newcommand{\gvec}{\boldsymbol{\mathsf{g}}} \newcommand{\zerovec}{\boldsymbol{\mathsf{0}}} \newcommand{\norm}[1]{\bigl\lVert{#1}\bigr\rVert} \newcommand{\transpose}[1]{{#1}^\top} \DeclareMathOperator{\rank}{rank} $
Refer to course Notebook 09.
'''Read a reaction mechanism and create data structures'''
from chen_3170.toolkit import reaction_mechanism # replace with your code
# build the stoichiometric matrix
(species, reactions, stoic_mtrx, dummy, dummy) = reaction_mechanism('data/water-gas-shift-rxn.txt', shuffle=False)
print(species)
from chen_3170.help import print_reactions
print_reactions(reactions)
['H2', 'H', 'CO2', 'OH', 'H2O', 'CO'] r0 : H2O + CO <=> CO2 + H2 r1 : H2O + H <=> H2 + OH r2 : OH + CO <=> CO2 + H n_reactions = 3
'''Check the stoichiometric matrix'''
from chen_3170.help import plot_matrix
plot_matrix(stoic_mtrx, title='Stoichiometric Matrix')
import numpy as np
np.set_printoptions(precision=3,threshold=100,edgeitems=5)
print('stoic_mtrx=\n',stoic_mtrx)
matrix shape = (3, 6)
stoic_mtrx= [[ 1. 0. 1. 0. -1. -1.] [ 1. -1. 0. 1. -1. 0.] [ 0. 1. 1. -1. 0. -1.]]
Refer to course Notebook 07.
'''Build the full-rank sub-mechanism reactions list'''
from chen_3170.toolkit import sub_mechanisms # replace with your code
sub_mechanisms = sub_mechanisms( species, reactions, stoic_mtrx )
# reactions = 3 # species = 6 rank of S = 2 # of all possible sub_mechanisms = 3 # of full-rank sub_mechanisms = 3
'''Top reaction sub-mechanism'''
sub_mechanism_1 = sub_mechanisms[0]
print(sub_mechanism_1)
[(0, 1), ['H2O + CO <=> CO2 + H2', 'H2O + H <=> H2 + OH'], array([[ 1., 0., 1., 0., -1., -1.], [ 1., -1., 0., 1., -1., 0.]]), 10.0]
'''Top reaction sub-mechanism stoichiometric matrix'''
stoic_mtrx_1 = sub_mechanism_1[2]
#( dummy, dummy, stoic_mtrx_1 ) = reaction_mechanism( reactions = reactions_1 ) # taking advantage of this function
#plot_matrix(stoic_mtrx_1, title='Sub-Mech 1')
print('S_1=\n',stoic_mtrx_1)
S_1= [[ 1. 0. 1. 0. -1. -1.] [ 1. -1. 0. 1. -1. 0.]]
'''Another top reaction sub-mechanism'''
sub_mechanism_2 = sub_mechanisms[1]
print(sub_mechanism_2)
[(0, 2), ['H2O + CO <=> CO2 + H2', 'OH + CO <=> CO2 + H'], array([[ 1., 0., 1., 0., -1., -1.], [ 0., 1., 1., -1., 0., -1.]]), 10.0]
'''Another top reaction sub-mechanism stoichiometric matrix'''
stoic_mtrx_2 = sub_mechanism_2[2]
#( dummy, dummy, stoic_mtrx_2 ) = reaction_mechanism( reactions = reactions_2 ) # taking advantage of this function
#plot_matrix(stoic_mtrx_1, title='Sub-Mech 2')
print('S_2=\n',stoic_mtrx_2)
S_2= [[ 1. 0. 1. 0. -1. -1.] [ 0. 1. 1. -1. 0. -1.]]
Refer to course Notebook 07.
'''Use given production rate'''
import numpy as np
# ['OH', 'CO', 'CO2', 'H2', 'H', 'H2O']
g_exp = dict()
g_exp['H'] = -1
g_exp['H2'] = 3
g_exp['OH'] = 1
g_exp['H2O'] = -3
g_exp['CO'] = -2
g_exp['CO2'] = 2
g_vec = np.array( [g_exp[spc] for spc in species] )
Here, let's compute $\rvec_1$ for $ \Smtrx_1\,\Smtrx_1^\top\,\rvec_1 = \Smtrx_1\,\gvec . $
'''Compute the LS reaction rates for random species production rates'''
# build A x = b LS problem
a_mtrx = stoic_mtrx_1 @ stoic_mtrx_1.transpose() # A = S ST, A is the normal matrix
b_vec = stoic_mtrx_1 @ g_vec # b = S g
from chen_3170.toolkit import lu_factorization # replace with your code
# matrix LU factorization of A, the normal matrix
(P,L,U,s_rank) = lu_factorization( a_mtrx, 'partial' ) # matrix is full rank; partial pivoting works
assert s_rank == np.linalg.matrix_rank(stoic_mtrx_1)
from chen_3170.help import forward_solve
from chen_3170.toolkit import backward_solve # replace with your code
# solve the LS problem: A x = b
y_vec = forward_solve( L, P @ b_vec) # L y = P b
x_vec = backward_solve( U, y_vec) # U x = y
assert np.linalg.norm(x_vec - np.linalg.solve(a_mtrx,b_vec)) < 1e-12
print('species production rates g_vec =',g_vec)
r_vec_1 = x_vec # r_1 = x
print('reaction rates r_vec_1=',r_vec_1)
residual_vec_1 = g_vec - stoic_mtrx_1.transpose() @ r_vec_1
print('residual norm g - ST_1 r_1 = %8.5e'%np.linalg.norm(residual_vec_1))
species production rates g_vec = [ 3 -1 2 1 -3 -2] reaction rates r_vec_1= [2. 1.] residual norm g - ST_1 r_1 = 0.00000e+00
'''Which reaction rates are these?'''
from chen_3170.help import print_reactions
print_reactions(sub_mechanism_1[1])
print(sub_mechanism_1[0])
print('')
print_reactions(reactions)
r0 : H2O + CO <=> CO2 + H2 r1 : H2O + H <=> H2 + OH n_reactions = 2 (0, 1) r0 : H2O + CO <=> CO2 + H2 r1 : H2O + H <=> H2 + OH r2 : OH + CO <=> CO2 + H n_reactions = 3
Here, let's compute $\rvec_2$ for $ \Smtrx_2\,\Smtrx_2^\top\,\rvec_2 = \Smtrx_2\,\gvec . $
'''Compute the LS reaction rates for random species production rates'''
# build A x = b LS problem
a_mtrx = stoic_mtrx_2 @ stoic_mtrx_2.transpose() # A = S ST, A is the normal matrix
b_vec = stoic_mtrx_2 @ g_vec # b = S g
from chen_3170.toolkit import lu_factorization # replace with your code
# matrix LU factorization of A
(P,L,U,s_rank) = lu_factorization( a_mtrx, 'partial' ) # matrix is full rank; partial pivoting works
assert s_rank == np.linalg.matrix_rank(stoic_mtrx_2)
from chen_3170.help import forward_solve
from chen_3170.toolkit import backward_solve # replace with your code
# solve the LS problem: A x = b
y_vec = forward_solve( L, P @ b_vec )
x_vec = backward_solve( U, y_vec )
assert np.linalg.norm(x_vec - np.linalg.solve(a_mtrx,b_vec)) < 1e-12
print('species production rates g_vec =',g_vec)
r_vec_2 = x_vec # r_2 = x
print('reaction rates r_vec_2=',r_vec_2)
residual_vec_2 = g_vec - stoic_mtrx_2.transpose() @ r_vec_2
print('residual norm g - ST_2 r_2 = %8.5e'%np.linalg.norm(residual_vec_2))
species production rates g_vec = [ 3 -1 2 1 -3 -2] reaction rates r_vec_2= [ 3. -1.] residual norm g - ST_2 r_2 = 0.00000e+00
'''Which reaction rates are these?'''
from chen_3170.help import print_reactions
print_reactions(sub_mechanism_2[1])
print(sub_mechanism_2[0])
print('')
print_reactions(reactions)
r0 : H2O + CO <=> CO2 + H2 r1 : OH + CO <=> CO2 + H n_reactions = 2 (0, 2) r0 : H2O + CO <=> CO2 + H2 r1 : H2O + H <=> H2 + OH r2 : OH + CO <=> CO2 + H n_reactions = 3
Note that the sub-mechanisms rates $\rvec_1$ and $\rvec_2$ are different, but the corresponding residuals $\gvec - \Smtrx_1^\top\,\rvec_1$, and $\gvec - \Smtrx_2^\top\,\rvec_2$, are equal!
print('r_vec_1 - r_vec_2 =',r_vec_1 - r_vec_2)
np.set_printoptions(precision=3)
print('g - ST_1 r_1 - (g - ST_2 r_2) =',(g_vec-stoic_mtrx_1.transpose()@r_vec_1)-(g_vec-stoic_mtrx_2.transpose()@r_vec_2))
r_vec_1 - r_vec_2 = [-1. 2.] g - ST_1 r_1 - (g - ST_2 r_2) = [0. 0. 0. 0. 0. 0.]
Let us verify the previous residual result for all sub-mechanims. Here we will compute $\Smtrx_k$ for $k=1,\ldots,$# of full-rank sub-mechanisms and the corresponding LS reaction rates $\rvec_k$ and residuals
\begin{equation*} \gvec - \Smtrx_k^\top\,\rvec_k . \end{equation*}'''Compute the LS residual for all full-rank reaction sub-mechanisms'''
import numpy as np
sub_mech_rxn_rates_mtrx = np.zeros( (s_rank,len(sub_mechanisms)) ) # column-wise storage of all r_k's
sub_mech_residuals_mtrx = np.zeros( (len(species),len(sub_mechanisms)) ) # column-wise storage of all g - ST_k r_k
np.set_printoptions(precision=2)
for sub_mech in sub_mechanisms:
stoic_mtrx_k = sub_mech[2] # view of stoichiometric matrix for sub-mechanism
# assemble LS problem matrix and right side
a_mtrx = stoic_mtrx_k @ stoic_mtrx_k.transpose() # A = S ST, A is the normal matrix
b_vec = stoic_mtrx_k @ g_vec # b = S g
# matrix LU factorization of A
(P,L,U,s_rank) = lu_factorization( a_mtrx, 'partial' ) # matrix is full rank; partial pivoting works
assert s_rank == np.linalg.matrix_rank(stoic_mtrx_k)
# solve the LS problem: A x = b
y_vec = forward_solve( L, P @ b_vec) # L y = P b
x_vec = backward_solve( U, y_vec) # U x = y
assert np.linalg.norm(x_vec - np.linalg.solve(a_mtrx,b_vec)) < 1e-12
r_vec = x_vec # r = x
sub_mech_rxn_rates_mtrx[:,sub_mechanisms.index(sub_mech)] = r_vec
print('rxns ids =',sub_mech[0])
print(' r =', r_vec)
rate_norm = np.linalg.norm(r_vec)
print(' ||r|| = %8.3e'%rate_norm)
residual = g_vec - stoic_mtrx_k.transpose() @ r_vec
sub_mech_residuals_mtrx[:,sub_mechanisms.index(sub_mech)] = residual
print('g - ST r =',residual)
print('||g - ST r|| = %8.3e'%np.linalg.norm(residual))
print('')
rxns ids = (0, 1) r = [2. 1.] ||r|| = 2.236e+00 g - ST r = [0. 0. 0. 0. 0. 0.] ||g - ST r|| = 0.000e+00 rxns ids = (0, 2) r = [ 3. -1.] ||r|| = 3.162e+00 g - ST r = [0. 0. 0. 0. 0. 0.] ||g - ST r|| = 0.000e+00 rxns ids = (1, 2) r = [3. 2.] ||r|| = 3.606e+00 g - ST r = [0. 0. 0. 0. 0. 0.] ||g - ST r|| = 0.000e+00
'''Plot the LS reaction rates for all full-rank sub-mechanisms'''
from matplotlib import pyplot as plt # import the pyplot function of the matplotlib package
fig, ax = plt.subplots(figsize=(20,6))
ax.plot(range(sub_mech_rxn_rates_mtrx.shape[1]), np.linalg.norm(sub_mech_rxn_rates_mtrx,axis=0),
'-.',color='black', marker='*',markersize=12)
plt.xticks(range(sub_mech_rxn_rates_mtrx.shape[1]),[sm[0] for sm in sub_mechanisms],rotation=60,fontsize=14)
ax.set_ylabel('Reaction Rate Norm',fontsize=16)
ax.set_xlabel('Reaction Sub-Mechanism',fontsize=16)
ax.xaxis.grid(True,linestyle='-',which='major',color='lightgrey',alpha=0.9)
fig.suptitle('Full-Rank (%i) LS Reaction Rates Magnitude of All Sub-Mechanisms'%(s_rank),fontsize=20)
plt.show()
'''Print the LS reaction rate with minimum norm (unique)'''
rate_min = np.linalg.norm(sub_mech_rxn_rates_mtrx,axis=0).min()
idx_min = np.argmin( np.linalg.norm(sub_mech_rxn_rates_mtrx,axis=0) )
print('sub-mechanism id. = %i'%idx_min)
print('minimum norm sub-mech rxn rate = %8.4e'%rate_min)
print('sub-mechanism rxn id =', sub_mechanisms[idx_min][0])
print_reactions(sub_mechanisms[idx_min][1])
print('sub-mechanism rxn score = %4.2f'%sub_mechanisms[idx_min][3])
sub-mechanism id. = 0 minimum norm sub-mech rxn rate = 2.2361e+00 sub-mechanism rxn id = (0, 1) r0 : H2O + CO <=> CO2 + H2 r1 : H2O + H <=> H2 + OH n_reactions = 2 sub-mechanism rxn score = 10.00
'''Plot the LS residual norms of all full-rank sub-mechanisms'''
from matplotlib import pyplot as plt # import the pyplot function of the matplotlib package
fig, ax = plt.subplots(figsize=(20,6))
ax.plot(range(sub_mech_residuals_mtrx.shape[1]), np.linalg.norm(sub_mech_residuals_mtrx,axis=0),
'-.',color='red',marker='*',markersize=12)
plt.xticks(range(sub_mech_residuals_mtrx.shape[1]),[smr[0] for smr in sub_mechanisms],rotation=60,fontsize=14)
ax.set_ylabel('Residual Norm',fontsize=16)
ax.set_xlabel('Reaction Sub-Mechanisms',fontsize=16)
ax.xaxis.grid(True,linestyle='-',which='major',color='lightgrey',alpha=0.9)
fig.suptitle('Full-Rank LS Residuals of All Sub-Mechanisms w/ Rank = '+str(s_rank),fontsize=20)
plt.show()
Despite all the foregoing development, we have not solved the original problem yet, namely
\begin{equation*} \Smtrx^\top\,\rvec = \gvec . \end{equation*}However this problem is related to the series of sub-mechanisms we just analyzed. To see this, assemble the matrix of sub-mechanism reaction rates mapped into the original reaction dimension $m$ by completing the additional entries with zeros. We denote this vector mapping $\rvec_k\longrightarrow\hat{\rvec}_k$ and the vectors $\hat{\rvec}_k$ as basic vectors.
Then we show that all basic vectors produce the same residual vectors, that is:
\begin{equation*} \gvec - \Smtrx^\top\,\hat{\rvec}_k \end{equation*}is invariant for any $\hat{\rvec}_k$ and has the same norm as the residuals $\gvec - \Smtrx_k^\top\,\rvec_k$. All residuals being equal, we choose the rank-deficient least-squares basic solution as the $\hat{\rvec}_k$ of minimum norm.
'''Map r_k --> \hat{r} rates'''
mech_rxn_rates_mtrx = np.zeros((len(reactions),len(sub_mechanisms)))
for sm in sub_mechanisms:
sm_idx = sub_mechanisms.index(sm)
rxn_idxs = sm[0]
mech_rxn_rates_mtrx[rxn_idxs,sm_idx] = sub_mech_rxn_rates_mtrx[:,sm_idx] # map
'''Show the mapping in the form of a matrix'''
from chen_3170.help import plot_matrix
plot_matrix(mech_rxn_rates_mtrx)
np.set_printoptions(precision=2,threshold=100,edgeitems=5)
print(mech_rxn_rates_mtrx)
matrix shape = (3, 3)
[[ 2. 3. 0.] [ 1. 0. 3.] [ 0. -1. 2.]]
Compute $\gvec - \Smtrx^\top\,\hat{\rvec}_k$ for all $k$.
'''Subtract above from g column by column'''
np.set_printoptions(precision=3)
for j in range(mech_rxn_rates_mtrx.shape[1]):
tmp = g_vec - (stoic_mtrx.transpose() @ mech_rxn_rates_mtrx)[:,j]
print('residual k = %i'%j,tmp)
residual k = 0 [0. 0. 0. 0. 0. 0.] residual k = 1 [0. 0. 0. 0. 0. 0.] residual k = 2 [0. 0. 0. 0. 0. 0.]
'''Plot of the input species production rate'''
from matplotlib import pyplot as plt # import the pyplot function of the matplotlib package
fig, ax = plt.subplots(figsize=(20,6))
ax.bar(range(len(species)), g_vec, color='blue')
plt.xticks(range(len(species)),species,rotation=60,fontsize=14)
ax.set_ylabel('Species Production Rate',fontsize=16)
ax.set_xlabel('Species',fontsize=16)
ax.xaxis.grid(True,linestyle='-',which='major',color='lightgrey',alpha=0.9)
fig.suptitle('Input Species Production Rates ('+str(len(species))+')',fontsize=20)
plt.show()
Select as basic solution of the rank-deficient LS problem:
\begin{equation*} \min\limits_k \norm{\hat{\rvec}_k} \quad\ \forall \quad\ k=1,\ldots,\text{# of full-rank sub-mechanisms} . \end{equation*}'''Basic least-squares reaction rates'''
from matplotlib import pyplot as plt # import the pyplot function of the matplotlib package
fig, ax = plt.subplots(figsize=(20,6))
ax.bar(range(len(reactions)), mech_rxn_rates_mtrx[:,idx_min], color='orange')
plt.xticks(range(len(reactions)),reactions,rotation=60,fontsize=14)
ax.set_ylabel('Reaction rate',fontsize=16)
ax.set_xlabel('Reaction',fontsize=16)
ax.xaxis.grid(True,linestyle='-',which='major',color='lightgrey',alpha=0.9)
fig.suptitle('Basic LS Reaction Rates (%s)'%str(s_rank),fontsize=20)
plt.show()
To compute the reaction rates vector $\rvec$ for a given species production vector $\gvec$ we need to solve:
\begin{equation*} \Smtrx^\top\,\rvec = \gvec . \end{equation*}Since $\Smtrx$ is rank deficient, a unique mininum-norm least squares solution, $\rvec_\text{MNLS}$, exists as follows,
\begin{equation*} \lim\limits_{\delta\rightarrow 0}\,\min\limits_\rvec \bigl( \norm{\gvec - \Smtrx^\top\,\rvec_\text{MNLS}}^2 + \delta\norm{\rvec_\text{MNLS}}^2 \bigr) \quad\ \forall \quad\ \rvec. \end{equation*}This solution is obtained by solving:
\begin{equation*} \lim\limits_{\delta\rightarrow 0}\,\bigl(\Smtrx\,\Smtrx^\top + \delta\Imtrx\bigr) \rvec_\text{MNLS} = \Smtrx\,\gvec , \end{equation*}where $\Smtrx\,\Smtrx^\top$ is square, symmetric and singular. The least-squares problem is just $\Amtrx\,\xvec=\bvec$ with $\Amtrx = \Smtrx\,\Smtrx^\top + \delta\Imtrx$ and $\bvec = \Smtrx\,\gvec$.
'''Compute the LS reaction rates for random species production rates'''
# Overdetermined
assert stoic_mtrx.transpose().shape[0] > stoic_mtrx.transpose().shape[1]
# build A x = b LS problem
delta = 1e-3
a_mtrx = stoic_mtrx @ stoic_mtrx.transpose() + delta * np.eye(3) # A = S ST + delta I, A is the normal matrix
b_vec = stoic_mtrx @ g_vec # b = S g
from chen_3170.toolkit import lu_factorization # replace with your code
# matrix LU factorization of A, the normal matrix
(P,L,U,a_rank) = lu_factorization( a_mtrx, 'partial' ) # matrix is full rank; partial pivoting works
assert a_rank == np.linalg.matrix_rank(stoic_mtrx@stoic_mtrx.transpose()+np.eye(3))
from chen_3170.help import forward_solve
from chen_3170.toolkit import backward_solve # replace with your code
# solve the LS problem: A x = b
y_vec = forward_solve( L, P @ b_vec) # L y = P b
x_vec = backward_solve( U, y_vec) # U x = y
assert np.linalg.norm(x_vec - np.linalg.solve(a_mtrx,b_vec)) < 1e-12
print('species production rates g_vec =',g_vec)
r_vec = x_vec # r = x
print('reaction rates r_vec=',r_vec)
rate_norm = np.linalg.norm(r_vec)
print(' ||r|| = %8.3e'%rate_norm)
residual_vec = g_vec - stoic_mtrx.transpose() @ r_vec
print('||g - ST r|| + ||r|| = %8.5e'%(np.linalg.norm(residual_vec)+rate_norm))
species production rates g_vec = [ 3 -1 2 1 -3 -2] reaction rates r_vec= [1.666 1.333 0.333] ||r|| = 2.160e+00 ||g - ST r|| + ||r|| = 2.16077e+00
'''Minimum norm least-squares reaction rates'''
from matplotlib import pyplot as plt # import the pyplot function of the matplotlib package
fig, ax = plt.subplots(figsize=(20,6))
ax.bar(range(len(reactions)), r_vec, color='orange')
plt.xticks(range(len(reactions)),reactions,rotation=60,fontsize=14)
ax.set_ylabel('Reaction rate',fontsize=16)
ax.set_xlabel('Reaction',fontsize=16)
ax.xaxis.grid(True,linestyle='-',which='major',color='lightgrey',alpha=0.9)
fig.suptitle('Minimum Norm LS Reaction Rates (%s)'%str(s_rank),fontsize=20)
plt.show()