ChEn-3170: Computational Methods in Chemical Engineering Fall 2018 UMass Lowell; Prof. V. F. de Almeida 22Oct2018
$ \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/nox-rxn.txt')
print(species)
from chen_3170.help import print_reactions
print_reactions(reactions)
['NO3', 'NO', 'N2O4', 'N2O5', 'O2', 'NO2'] r0 : 2 N2O5 <=> 2 N2O4 + 1 O2 r1 : 1 N2O5 <=> 1 NO2 + 1 NO3 r2 : 1 NO2 + 1 NO3 <=> 1 N2O5 r3 : 1 NO3 <=> 1 O2 + 1 NO r4 : 1 NO + 1 N2O5 <=> 3 NO2 r5 : 2 NO2 <=> 1 N2O4 n_reactions = 6
'''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 = (6, 6)
stoic_mtrx= [[ 0. 0. 2. -2. 1. 0.] [ 1. 0. 0. -1. 0. 1.] [-1. 0. 0. 1. 0. -1.] [-1. 1. 0. 0. 1. 0.] [ 0. -1. 0. -1. 0. 3.] [ 0. 0. 1. 0. 0. -2.]]
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 = 6 # species = 6 rank of S = 4 # of all possible sub_mechanisms = 15 # of full-rank sub_mechanisms = 9
'''Top reaction sub-mechanism'''
sub_mechanism_1 = sub_mechanisms[0]
print(sub_mechanism_1)
[(0, 3, 4, 5), ['2 N2O5 <=> 2 N2O4 + 1 O2', '1 NO3 <=> 1 O2 + 1 NO', '1 NO + 1 N2O5 <=> 3 NO2', '2 NO2 <=> 1 N2O4'], array([[ 0., 0., 2., -2., 1., 0.], [-1., 1., 0., 0., 1., 0.], [ 0., -1., 0., -1., 0., 3.], [ 0., 0., 1., 0., 0., -2.]]), 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= [[ 0. 0. 2. -2. 1. 0.] [-1. 1. 0. 0. 1. 0.] [ 0. -1. 0. -1. 0. 3.] [ 0. 0. 1. 0. 0. -2.]]
'''Another top reaction sub-mechanism'''
sub_mechanism_2 = sub_mechanisms[1]
print(sub_mechanism_2)
[(0, 1, 3, 4), ['2 N2O5 <=> 2 N2O4 + 1 O2', '1 N2O5 <=> 1 NO2 + 1 NO3', '1 NO3 <=> 1 O2 + 1 NO', '1 NO + 1 N2O5 <=> 3 NO2'], array([[ 0., 0., 2., -2., 1., 0.], [ 1., 0., 0., -1., 0., 1.], [-1., 1., 0., 0., 1., 0.], [ 0., -1., 0., -1., 0., 3.]]), 8.928571428571429]
'''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= [[ 0. 0. 2. -2. 1. 0.] [ 1. 0. 0. -1. 0. 1.] [-1. 1. 0. 0. 1. 0.] [ 0. -1. 0. -1. 0. 3.]]
Refer to course Notebook 07.
'''Assume a species production rate as random'''
import numpy as np
g_vec = np.random.random(len(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 = [0.51 0.365 0.306 0.67 0.468 0.915] reaction rates r_vec_1= [ 0.183 -0.038 -0.253 -0.682] residual norm g - ST_1 r_1 = 1.20229e+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 : 2 N2O5 <=> 2 N2O4 + 1 O2 r1 : 1 NO3 <=> 1 O2 + 1 NO r2 : 1 NO + 1 N2O5 <=> 3 NO2 r3 : 2 NO2 <=> 1 N2O4 n_reactions = 4 (0, 3, 4, 5) r0 : 2 N2O5 <=> 2 N2O4 + 1 O2 r1 : 1 N2O5 <=> 1 NO2 + 1 NO3 r2 : 1 NO2 + 1 NO3 <=> 1 N2O5 r3 : 1 NO3 <=> 1 O2 + 1 NO r4 : 1 NO + 1 N2O5 <=> 3 NO2 r5 : 2 NO2 <=> 1 N2O4 n_reactions = 6
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 = [0.51 0.365 0.306 0.67 0.468 0.915] reaction rates r_vec_2= [-0.158 0.341 0.303 0.088] residual norm g - ST_2 r_2 = 1.20229e+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 : 2 N2O5 <=> 2 N2O4 + 1 O2 r1 : 1 N2O5 <=> 1 NO2 + 1 NO3 r2 : 1 NO3 <=> 1 O2 + 1 NO r3 : 1 NO + 1 N2O5 <=> 3 NO2 n_reactions = 4 (0, 1, 3, 4) r0 : 2 N2O5 <=> 2 N2O4 + 1 O2 r1 : 1 N2O5 <=> 1 NO2 + 1 NO3 r2 : 1 NO2 + 1 NO3 <=> 1 N2O5 r3 : 1 NO3 <=> 1 O2 + 1 NO r4 : 1 NO + 1 N2O5 <=> 3 NO2 r5 : 2 NO2 <=> 1 N2O4 n_reactions = 6
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 = [ 0.341 -0.379 -0.557 -0.77 ] g - ST_1 r_1 - (g - ST_2 r_2) = [-2.220e-16 -1.665e-16 0.000e+00 1.110e-16 0.000e+00 0.000e+00]
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, 3, 4, 5) r = [ 0.18 -0.04 -0.25 -0.68] ||r|| = 7.515e-01 g - ST r = [0.47 0.15 0.62 0.78 0.32 0.31] ||g - ST r|| = 1.202e+00 rxns ids = (0, 1, 3, 4) r = [-0.16 0.34 0.3 0.09] ||r|| = 4.910e-01 g - ST r = [0.47 0.15 0.62 0.78 0.32 0.31] ||g - ST r|| = 1.202e+00 rxns ids = (0, 1, 3, 5) r = [-0.07 0.25 0.22 -0.18] ||r|| = 3.827e-01 g - ST r = [0.47 0.15 0.62 0.78 0.32 0.31] ||g - ST r|| = 1.202e+00 rxns ids = (0, 1, 4, 5) r = [ 0.15 0.04 -0.22 -0.61] ||r|| = 6.611e-01 g - ST r = [0.47 0.15 0.62 0.78 0.32 0.31] ||g - ST r|| = 1.202e+00 rxns ids = (0, 2, 3, 4) r = [-0.16 -0.34 0.3 0.09] ||r|| = 4.910e-01 g - ST r = [0.47 0.15 0.62 0.78 0.32 0.31] ||g - ST r|| = 1.202e+00 rxns ids = (0, 2, 3, 5) r = [-0.07 -0.25 0.22 -0.18] ||r|| = 3.827e-01 g - ST r = [0.47 0.15 0.62 0.78 0.32 0.31] ||g - ST r|| = 1.202e+00 rxns ids = (0, 2, 4, 5) r = [ 0.15 -0.04 -0.22 -0.61] ||r|| = 6.611e-01 g - ST r = [0.47 0.15 0.62 0.78 0.32 0.31] ||g - ST r|| = 1.202e+00 rxns ids = (1, 3, 4, 5) r = [ 0.18 0.15 -0.07 -0.32] ||r|| = 3.994e-01 g - ST r = [0.47 0.15 0.62 0.78 0.32 0.31] ||g - ST r|| = 1.202e+00 rxns ids = (2, 3, 4, 5) r = [-0.18 0.15 -0.07 -0.32] ||r|| = 3.994e-01 g - ST r = [0.47 0.15 0.62 0.78 0.32 0.31] ||g - ST r|| = 1.202e+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. = 2 minimum norm sub-mech rxn rate = 3.8266e-01 sub-mechanism rxn id = (0, 1, 3, 5) r0 : 2 N2O5 <=> 2 N2O4 + 1 O2 r1 : 1 N2O5 <=> 1 NO2 + 1 NO3 r2 : 1 NO3 <=> 1 O2 + 1 NO r3 : 2 NO2 <=> 1 N2O4 n_reactions = 4 sub-mechanism rxn score = 8.93
'''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 equal to 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 = (6, 9)
[[ 0.18 -0.16 -0.07 0.15 -0.16 -0.07 0.15 0. 0. ] [ 0. 0.34 0.25 0.04 0. 0. 0. 0.18 0. ] [ 0. 0. 0. 0. -0.34 -0.25 -0.04 0. -0.18] [-0.04 0.3 0.22 0. 0.3 0.22 0. 0.15 0.15] [-0.25 0.09 0. -0.22 0.09 0. -0.22 -0.07 -0.07] [-0.68 0. -0.18 -0.61 0. -0.18 -0.61 -0.32 -0.32]]
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.472 0.149 0.622 0.783 0.323 0.311] residual k = 1 [0.472 0.149 0.622 0.783 0.323 0.311] residual k = 2 [0.472 0.149 0.622 0.783 0.323 0.311] residual k = 3 [0.472 0.149 0.622 0.783 0.323 0.311] residual k = 4 [0.472 0.149 0.622 0.783 0.323 0.311] residual k = 5 [0.472 0.149 0.622 0.783 0.323 0.311] residual k = 6 [0.472 0.149 0.622 0.783 0.323 0.311] residual k = 7 [0.472 0.149 0.622 0.783 0.323 0.311] residual k = 8 [0.472 0.149 0.622 0.783 0.323 0.311]
'''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()