Assume that many blood vessels in a region of the body can ...



ESTIMATION OF DYNAMIC KINETIC RATE CONSTANTSJeff Marsh1, Andreas A. Linninger21jmarsh29@uic.edu , 2linninge@uic.edu Department of Bioengineering, University of Illinois at ChicagoThis report is produced under the supervision of BIOE310 instructor Prof. Linninger.AbstractA system of three reacting species with known reaction stoichiometry is considered in this report. It is of interest to investigate the ability of Parameter estimation techniques for inferring the reaction rate values from experimental data. To provide validation with a measurable amount and frequency of experimental error, experimental results were artificially generated with a forward simulation. To avoid an optimization problem with ordinary differential equation constraints, we discretize the species balances. When discretizing the differential equations representing the species balances for each reactant, the prediction of the kinetic rates becomes a linear algebraic problem representing the inversion of the pseudoinverse of the experimental data. Reaction rate estimation will also be evaluated when subject to an applied error percent, in order to determine the reliability of estimating parameters subject to experimental noise. Keywords: Eigenvalue decomposition, reaction rate, kinetic reaction, parameter estimation1. IntroductionThis project is interested in investigating parameter estimation techniques for a dynamic kinetic reaction system, as described by Froment and Bischoff 1. The model is useful in scenarios where the scientist/engineer has a system with a series of coupled reactions, and where the concentrations at time-points are all which can be attained. With gaining insight into how the concentrations of the chemical species of interest evolve with the progression of time, it is possible to deduce the reaction rate constants which govern the conversion of species to each other. This estimation of reaction parameters is important for many applications. For example: when analyzing a colony of cells, where concentration of reactants and products can be measure in time, but not the reaction rates governing the conversion of reactants to products.2. MethodsThis report will focus on a representative dynamic kinetic reacting system of three chemical species, which start at initial values and react toward equilibrium; Figure 1 found in the Appendix section illustrate this design.2.1 Species EquationsA batch reactor of finite volume contains three species (A, B, C), which have initial concentrations (CA0, CB0, CC0) at t=0. The batch reactor is assumed to be well mixed and contain equivalent species concentrations throughout its volume, and also be uninfluenced by thermal fluctuations. The generation or consumption of any reacting species is time dependent and governed by the six reaction rates (kAB, kba, kBC, kcb, kCA, kac), as well as the concentration of all species at that time. At t > 0 the species begin to react, and the time dependent equations in component form can be denoted as:?CA?t=-kAB+kac*CA+kba*CB+kCA*CC (1)?CB?t=kAB*CA-kba+kBC*CB+kcb*CC (2)?CC?t=kac*CA+kBC*CB-kcb+kCA*CC (3)It is then possible to rewrite the component form equations (1), (2), and (3) in matrix form:-12065011430?CA?t?CB?t?CC?t=-kAB+kackbakCAkAB-kba+kBCkcbkackBC-kcb+kCACACBCC (4)00?CA?t?CB?t?CC?t=-kAB+kackbakCAkAB-kba+kBCkcbkackBC-kcb+kCACACBCC (4)When the reaction rates are known artificial measurements can be obtained for all reacting species by performing Eigenvalue Decomposition to the matrix proposed. To estimate reaction rates from experiments, an experimental dataset of changing species concentrations over a time interval was generated by forward simulation with known reaction rates. For the purposes of this report, synthetic experimental datasets can be generated by using simulations of concentration profiles from known reaction rate constants.2.2 Synthesis of Measurement Data bySimulationExperimental Data can be obtained from the reacting system within the batch reactor by taking samples at a set time increment. The batch reactor is well mixed, and therefore each sample taken will contain concentration information for each species at that sampled time. It is possible to use the concentration measurements for each species over the sampled time interval to estimate the reaction rate constants, which govern the reacting species concentration changes. To begin estimation of kinetic parameters, the parameter estimation method defined by Froment and Bischoff is applied to discretize the differential equations into algebraic difference equations1. Therefore the component equations (1), (2), and (3) can be integrated between t(i) and t(i+1), which yields algebraic equations (4), (5), and (6).-23685550380titi+1?CA=titi+1-kAB+kac?CA+kba?CB+kCA?CC?t (4)titi+1?CB=titi+1{ kAB?CA-kba+kBC?CB+kcb?CC}?t (5)titi+1?CC=titi+1 {kac?CA+kBC?CB-kcb+kCA?CC}?t (6)00titi+1?CA=titi+1-kAB+kac?CA+kba?CB+kCA?CC?t (4)titi+1?CB=titi+1{ kAB?CA-kba+kBC?CB+kcb?CC}?t (5)titi+1?CC=titi+1 {kac?CA+kBC?CB-kcb+kCA?CC}?t (6)After evaluating the integrals of each component equation respectively we obtain algebraic difference equations (7), (8), and (9). These equations can then be placed into matrix form by dividing out a column vector of the unknown reaction rates; providing a system of linear equations, (10), where all concentrations will be known from experimental data.-1762667716CAti+1-CAti=?t2[{-kAB+kac?CAti+1+kba?CBti+1+kCA?CCti+1} + {-kAB+kac?CAti+kba?CBti+kCA?CCti}] (7)CBti+1-CBti=?t2[{ kAB?CAti+1-kba+kBC?CBti+1+kcb?CCti+1} + { kAB?CAti-kba+kBC?CBti+kcb?CCti}] (8)CCti+1-CCti=?t2[{kac?CAti+1+kBC?CBti+1-kcb+kCA?CCti+1}+ {kac?CAti+kBC?CBti-kcb+kCA?CCti}] (9)CAti+1-CAti=?t2[{-kAB+kac?CAti+1+kba?CBti+1+kCA?CCti+1} + {-kAB+kac?CAti+kba?CBti+kCA?CCti}] (7)CBti+1-CBti=?t2[{ kAB?CAti+1-kba+kBC?CBti+1+kcb?CCti+1} + { kAB?CAti-kba+kBC?CBti+kcb?CCti}] (8)CCti+1-CCti=?t2[{kac?CAti+1+kBC?CBti+1-kcb+kCA?CCti+1}+ {kac?CAti+kBC?CBti-kcb+kCA?CCti}] (9)-17583539370?t2-CAti+1-CAti-CAti+1-CAti0CAti+1+CAti0-CBti+1-CBti0CAti+1+CAtiCBti+1+CBti CBti+1+CBtiCCti+1+CCti0-CBti+1-CBti0CCti+1+CCti0-CCti+1-CCti-CCti+1-CCtikABkackBCkbakCAkcb 10 ==CAti+1-CAtiCAti+1-CAtiCAti+1-CAti 00?t2-CAti+1-CAti-CAti+1-CAti0CAti+1+CAti0-CBti+1-CBti0CAti+1+CAtiCBti+1+CBti CBti+1+CBtiCCti+1+CCti0-CBti+1-CBti0CCti+1+CCti0-CCti+1-CCti-CCti+1-CCtikABkackBCkbakCAkcb 10 ==CAti+1-CAtiCAti+1-CAtiCAti+1-CAti 2.3 Artificially Added Measurement ErrorThe linear sum equations in the matrix form of equation (10) are all dependent on species concentration values between two time points under evaluation. Note that the linear set of equations (10) relates known rates to experimental data in the linear form of equation (11). The data measurements within the matrix represent the measured species concentrations at sampled times. The time points are presumed to be dictated by the experimental sampling rate, which is limited by the laboratory measuring equipment. Another limitation of laboratory measuring equipment is introduction of white noise into the data being measured. Examination of the effects of experimental error is also available for evaluation in this computer based experiment.Using the known kinetic rate constants to generate the species concentration profiles, which we can evaluate at any time points the sampling rate determines, it is then possible to add synthetic error to the experimental time points. Generation of synthetic error in MATLAB can be done by using a random number generator, and applying an upper and lower bound to obtain a specified range of error. The artificially generated errors can be added to the synthetically obtained species concentration data points, which then results in noisy data.The effect of the white noise in the experimental data on estimated reaction rates can be observed by implementing different amounts of error. In this study we will investigate the repercussions of adding 5%, 10%, and 20% error into our experimental data to observe it effect on the kinetic parameter analysis. This will also be demonstrative of faulty laboratory measuring equipment and how experimenters can be misled by inaccurate findings.Ak=C (11)2.4 Kinetic Parameter estimationAccording to equation (11) the difference of species concentrations at two distinct time points is equal to the product of the matrix containing summation equations of experimental species concentration data, and the unknown kinetic rate constant vector. Within each element of matrix A in equation (10) resides the sum of concentrations, where the concentration of a species at the first time point is added to the next; that is, at t(i) and t(i+1). All successive time points are separated by the same ?t, and the value is of interest because it may play an important role in the accuracy of the estimation of the unknown reaction rates. It can be noted that the matrix (10) which is generated by equations (7), (8), and (9) is underdetermined between t(i) and t(i+1). However, by evaluating the equations (7), (8), and (9) between t(i) and t(i+1), and allowing them to be the first 3 rows of the matrix, we can then reevaluate equations (7), (8), and (9) between t(i+1) and t(i+2), and set them as rows 4, 5, and 6. If we continue this process for many time intervals, we then have an overdetermined matrix. The overdetermined matrix now converts the problem into a pseudoinverse problem, where the unknown kinetic rates can be estimated with equation (12).k=(ATA)-1ATC (12)It can be seen that the k vector solution of equation (12) is the vector of supposed unknown reaction rates. It is noted that this explained method was intended to estimate the original defined reaction rates, from which the synthetic data of species concentrations at various times was generated.Therefore the accuracy of the method can be determined by resolving the original problem with the estimated k values, and plot the original concentration profiles against those of the estimated values. The effect of step size and time interval on the values may also be evaluated.3. ResultsIn order to evaluate the process described in the methods section for estimating unknown reaction rates in a three species reacting system within a constant volume batch reactor, initial values such as: initial concentrations(CA0, CB0, CC0), target kinetic rate constants(kAB, kBC, kba, kcb, kCA, kac), and a sampling rate(?t) must be chosen. To ensure adequate results will be obtained, six different experimental runs with differing initial reactor concentrations were done with three different sampling rates; noted in Table 2. The MATLAB code responsible for executing the procedure, Parameter_Estimation, is found in the MATLAB codes section of the appendix. This program is designed to generate synthetic experimental data for species concentrations using defined kinetic rate constants, and then estimate the target kinetic rate constants using the experimental data. The final aim is to graphically compare the concentration profiles obtained with the estimated kinetic rate constants to the experimentally sampled profiles.3.1 Experimental Runs Without ErrorIn order to apply the methods of parameter estimation we define the target kinetic rate constant values to be: kAB=4, kBC=3, kba=6.66667, kcb=4.5, kCA=4, kac=1.6; and initial species concentrations within the batch reactor to be pure A (CA0 =1, CB0=0 , CC0 =0). After defining initial states we rely on the function ConcProf, found in the MATLAB codes section of the appendix, to determine the concentration profiles CA(t), CB(t), and CC(t). To visualize the behavior of the concentration profiles CA(t), CB(t), and CC(t) with progression of time, we can observe a graphical plot of the original concentrations in Figure 2 of the appendix.The concentration profiles CA(t), CB(t), and CC(t) can now be used as experimentally obtained concentration data, which can be utilized to solve for the target kinetic rate constants. The concentration profiles must then be sampled to simulate experimental measurement data from the batch reactor. A sampling rate of ?t = 2 minutes/sample was used over the time interval t0 = 0 to tf = 15 minutes. The function EstRxnRates found in the MATLAB codes section of the appendix uses the defined sampling rate and time interval to perform the parameter estimation procedure on the experimentally obtained concentrations. The procedure mainly involves: discretization over defined time points and psuedoinverse evaluation.The function returns estimated kinetic rate constants of: kAB=2.59, kBC=1.84, kba=4.33, kcb=2.75, kCA=3.48, kac=1.39.The accuracy of the estimated kinetic rate constants can be evaluated by generating concentration profiles using the estimated kinetic rate constants, and then plot them with the original profiles. To implement this strategy the function ConcProfEst, found in the MATLAB codes section of the appendix, can be used. This function generates estimated concentration profiles through eigenvalue decomposition. Figure 4 below shows the resultant graphical representation of the original profiles and the estimated profiles. Upon observation of Figure 4, it can be noted that at the chosen ?t, there exists error in the estimated reaction rate constants even though the estimated concentration profiles were approximately accurate. The kinetic reaction rates obtained can be noted in Table 3 below.023938Figure 4. Graphical representation of the estimated parameters concentrations (CA(t)estimate, CB(t)estimate, and CC(t)estimate)with the target concentration behavior (CA(t), CB(t), and CC(t)), sampling rate of 1sample/2 minutes.0Figure 4. Graphical representation of the estimated parameters concentrations (CA(t)estimate, CB(t)estimate, and CC(t)estimate)with the target concentration behavior (CA(t), CB(t), and CC(t)), sampling rate of 1sample/2 minutes.092590Table 3. A table representing initial species concentrations of C0 = [1;0;0]. Target Kinetic reaction rates and estimated kinetic reaction rates can be observed.00Table 3. A table representing initial species concentrations of C0 = [1;0;0]. Target Kinetic reaction rates and estimated kinetic reaction rates can be observed.kiOriginalEstimatedkAB0.30.2920kac0.160.156kBC0.10.099kba0.10.097kCA0.10.097kcb0.150.1493.2 Runs with Variable Sample TimeAltering the chosen sampling rate ?t can have a significant impact on the estimation of kinetic rate constants. Figure 3 and Figure 5 are graphical representations of the same parameter estimation process, however instead ?t = 1 and ?t = 4 minutes/sample respectively. It can be noted that Figure 5 has the most amount of error, and Figure 3 has the least. Furthermore Table 2 in the appendix shows the impact of the different sampling rates on the estimated kinetic rate constants. Table 2 was obtained by recording the estimated kinetic parameters obtained by successively defining a new ?t and re-calling EstRxnRates. It is clear that the higher the frequency that samples are taken at, the less uncertainty exists in the measures, and accuracy in the calculations increases dramatically. It follows that the uncertainty in the data due to number of samples taken is proportional to the amount of error which will appear in resulting calculations.3.3 Runs Introducing Experimental ErrorIntroducing the artificial error into the experimental data of species concentration within the batch reactor at different time points can also have a profound effect on the estimation of kinetic rate constants. In this report we analyzed the effects of parameter estimation with induced 5%, 10%, and 20% artificial error added into the experimental species concentration measurements. The error was generated randomly between specified upper and lower bounds denoted, +/- 5% for example, then added into the data at every measurement. The function called to do this was Artificial_error, and the outcomes for the estimated reaction rates for the error data was obtained by EstRxnRatesErr, both can be found in the MATLAB codes section of the appendix. The resultant kinetic rate constants were placed into Table 1. Graphical representations Figures 7-9 show the plotted concentration data obtained using the error kinetic rate constants. As the values for the estimated kinetic rate constants show in Table 1, or the behavior noted in Figures 7-9, there is a gradual increase in the error for the kinetic rate constants obtained from data with increasing measurement error percentage. Below Figure 7 shows the graphical representation of 5% error plot with true concentration curves, and concentration curves generated from error influenced kinetic rates. The 5% measurement error effect on the estimated kinetic rates can be noted in Table 4.5729459750Figure 7. This figure shows the synthetic experimental data for concentration species with no error, then with induced artificial error of 5%, then the concentration profiles obtained with the estimated reaction rates of the error data.0Figure 7. This figure shows the synthetic experimental data for concentration species with no error, then with induced artificial error of 5%, then the concentration profiles obtained with the estimated reaction rates of the error data.1339854960Table 4. A table representing initial species concentrations of C0 = [1;0;0]. Target Kinetic reaction rates and estimated kinetic reaction rates influenced by 5% measurement error can be observed.00Table 4. A table representing initial species concentrations of C0 = [1;0;0]. Target Kinetic reaction rates and estimated kinetic reaction rates influenced by 5% measurement error can be observed.kiOriginalEstimated 5% ErrorkAB0.30.2797kac0.160.1565kBC0.10.0117kba0.10.1431kCA0.10.0127kcb0.150.09334. Discussion4.1 Parameter Estimation and Initial ReactorConcentrationsAlthough in the results section we focused on exploring the effectiveness of estimating kinetic rate constants for one concentration, with an initial concentration of pure species A; we can further explore the technique by applying it to multiple initial batch reactor concentrations. Although initial species concentrations within the batch reactor can be altered, the kinetic rate constants determine that the species will reach the same equilibrium. To observe the evolution toward equilibrium of different concentration profiles, we can plot them in a reaction triangle. The reaction triangle notes the state space of two concentrations, chosen here to be CA(t) and CB(t). To generate an error free reaction triangle, we can use the MATLAB function PlotRxnTri found in the MATLAB Codes section of the appendix. Figure 6 displays the state space reaction triangle, starting with initial batch reactor species contents from pure concentrations of each of the species. Note that because each initial reactor concentration progresses toward equilibrium, the estimated kinetic rate constants must be nearly identical. Table 2 in the appendix shows the effectiveness of the parameter estimation process for six different experimental runs, each starting with a unique initial reactor concentration. The consistency of the program executing the parameter estimation process is clear, as the estimated kinetic rate constants are identical for each experimental run.4.2 Parameter estimation and Sampling RateInvestigating the effect of the experimental sampling rate of the species concentrations is of interest to the process. If there is an optimal sampling rate to achieve nearly exact estimations of the kinetic rate constants, then we should be able to discover it by testing different sampling rates. Figure 3, Figure 4, and Figure 5 in the appendix show that by increasing the sampling rate, the kinetic rate constant estimates become increasingly more exact to the true values. This corresponds to the amount of measurement error that exists in each sample taken. For this reason, the more samples taken, the easier it is to eliminate the error during the parameter estimation process. Conversely when the sampling rate is to slow, the measurement error among the samples will significantly disrupt the accuracy of the reaction rate estimates. Table 2 also denotes the same conclusive effect of the sampling rate size on the estimated kinetic rate constants. The table also shows that the sampling rate has pronounced effect on the estimation of the kinetic rate constants, while initial concentrations have almost none. It is ultimately desirable from this conclusion therefore, to have a rapid sampling rate; however with modern technology it is not always possible to achieve this goal.4.3 Simulations Incorporating ArtificialErrorThis process can also be simulated by adding experimental error to the synthetic data of the species concentrations. To investigate the effect of measurement error, MATLAB functions Artificial_error and EstRxnRates were called upon in the main script and are found in the MATLAB Codes section of the appendix. Figure 7, Figure 8, and Figure 9 are graphical plots of concentration profiles which show experimental runs where the effect of 5%, 10%, and 20% added experimental error respectively. Using the MATLAB function PlotRxnTriERR found in the MATLAB Codes section of the appendix, we generated Figure 10, Figure 11, and Figure 12 to show the reaction triangles for the added experimental errors to the concentration profiles progression to equilibrium for 5%, 10%, and 20% error respectively. The resultant influence of the measurement error on the parameter estimation process is that as the experimental error increases, the accuracy of the kinetic rate constants decreases. To further support this claim we can observe Table 1 in the appendix section. Table 1 contains the results obtained from MATLAB in a series of experimental runs, where the estimated rate constants were calculated with artificial measurement error of 5%, 10%, and 20%. The information contained there supports what is visualized in Figures 7-12, which expresses the trend that as the amount of error in the measurements increases, the accuracy of the estimated kinetic parameters will decrease.4.4 Lagrangian InterpolationOccasionally concentration measurement instruments only allow for a small number of data measures to be taken during the course of the reaction. In cases such as this it is difficult to apply the kinetic parameter estimation process and achieve high accuracy. Another application which can be used to interpret the data is Lagrange Interpolation. In this process, Lagrange polynomials are used to fit splines between the points of the original experimental data. The splines between points allow for estimation of a polynomial function which would fit them together, relative to all data points. In order to evaluate the method of Lagrangian interpolation, we assume our measurement instruments limit us to only recording four data points. A code is generated to utilize proper Lagrangian polynomials for four data points can be noted as LagrangeInt_4pt in the MATLAB codes portion of the appendix. Figure 13 is a graphical representation of the resulting splines generated by this method. The resultant lines of the species concentration changes are reminiscent of the true curves, however still illicit a variable amount of error. The kinetic rate constants obtained from analysis of the Lagrangian interpolated data can be observed in Table 5. These values of the rate constants also have a corresponding amount of error proportional to the inaccuracy of the concentration data curves. 02389038Figure 13. This plot shows the effectiveness of Lagrangian Interpolation to recreate concentration profiles from only 4 measured data points. The estimated profile using this supertemporal estimation method proves to be an adequate estimator. 0Figure 13. This plot shows the effectiveness of Lagrangian Interpolation to recreate concentration profiles from only 4 measured data points. The estimated profile using this supertemporal estimation method proves to be an adequate estimator. ?-111760-77147Table 5. A table representing initial species concentrations of C0 = [1;0;0]. Target Kinetic reaction rates and estimated kinetic reaction rates influenced by 5% measurement error can be observed.00Table 5. A table representing initial species concentrations of C0 = [1;0;0]. Target Kinetic reaction rates and estimated kinetic reaction rates influenced by 5% measurement error can be observed.kiOriginalEstimated LagrangekAB0.30.1973kac0.160.1128kBC0.10.7033kba0.1-1.7565kCA0.13.0197kcb0.15-1.81995. ConclusionThe method for estimating reaction rates is valid process which can be used to accurately deduce unknown reaction rates from experimental data of species concentration. Otherwise, as demonstrated in the results section, this method can allow for simulations of reacting systems by generating synthetic data and using it to obtain reaction rates, without ever having to run real chemical experiments. When collecting data during a chemical experiment, the sampling rate chosen can have a tremendous influence on the estimated kinetic parameters. White noise can be added to synthetic experimental data to emulate experimental error encountered with recording instruments in the laboratory; the estimated rate constants are affected consequentially.Intellectual PropertyBiological and physiological data and some modeling procedures provided to you from Dr. Linninger’s lab are subject to IRB review procedures and Intellectual property procedures.Therefore, the use of these data and procedures are limited to the coursework only. Publications need to be approved and require joint authorship with staff of Dr. Linninger’s lab.6. ReferencesG. Froment, K Bischoff. Chemical Reaction Analysis and Design, Wiley Series & Sons, 1990.Andreas A Linninger (2014). Dynamic Kinetic Reaction, Eigenvalue decomposition, pseudoinverse analysis,lecture note.S. Kim. Dynamic Kinetic reaction, Department of Bioengineering, University of Illinois at Chicago.6. AppendixA. Experimental Design-6711586865Figure 1. A model of the proposed experiment: a batch reactor with three chemical species (A, B, C) at initial concentrations (CA0, CB0, CC0) react governed by the reaction rates shown (kAB kac kBC kba kCA kcb).0Figure 1. A model of the proposed experiment: a batch reactor with three chemical species (A, B, C) at initial concentrations (CA0, CB0, CC0) react governed by the reaction rates shown (kAB kac kBC kba kCA kcb).B. Parameter EstimationFigure 2 is a graphical representation of the behavior of species A, B, and C as they react from initial concentrations to equilibrium within the batch reactor.-431802147570Figure 2. This figure is a graphical representation of the three species A, B, and C reacting to equilibrium from the initial batch reactor concentration of pure species A.0Figure 2. This figure is a graphical representation of the three species A, B, and C reacting to equilibrium from the initial batch reactor concentration of pure species A.B. Concentration Profile Plots and time-stepFigure 3, Figure 4, and Figure 5 are graphical plots containing the actual behavior of the reacting species in time to equilibrium, which are compared to concentration profiles generated from the estimated kinetic parameters. The estimated kinetic rates are evaluated at different sampling rates; ?t = 1, 2, and 4 minutes per sample. There is a clear relationship between high accuracy with a lower ?t between experimental measures, as noted in Figure 3 and Figure 4, and less accuracy with a larger time between experimental measures, as noted in Figure 5.6038517013Figure 3. Graphical representation of the estimated parameters concentrations (CA(t)estimate, CB(t)estimate, and CC(t)estimate)with the target concentration behavior (CA(t), CB(t), and CC(t)), sampling rate 1 sample/minute.0Figure 3. Graphical representation of the estimated parameters concentrations (CA(t)estimate, CB(t)estimate, and CC(t)estimate)with the target concentration behavior (CA(t), CB(t), and CC(t)), sampling rate 1 sample/minute.-2540-6985Figure 5. Graphical representation of the estimated parameters concentrations (CA(t)estimate, CB(t)estimate, and CC(t)estimate)with the target concentration behavior (CA(t), CB(t), and CC(t)), sampling rate of 1sample/4 minutes.Figure 5. Graphical representation of the estimated parameters concentrations (CA(t)estimate, CB(t)estimate, and CC(t)estimate)with the target concentration behavior (CA(t), CB(t), and CC(t)), sampling rate of 1sample/4 minutes.C. Reaction Triangle for Estimated ConcentrationsFor the pure experimental data, it is possible to gain insight from the reaction triangle of species concentrations. Figure 6 represents species A plotted against species B, with the hypotenuse of the reaction triangle being the boundary, and can be thought of as pure CC(t). As the legend shows, a pure concentration of CA(t), CB(t), and CC(t) were used as three distinct initial concentrations, and ?t=1 minutes/sample. All profiles converge to the same equilibrium regardless of initial concentrations.02256288Figure 6. Reaction triangle showing the advancement toward the equilibrium point for initial batch reactor concentrations of: pure A, pure B, and pure C.0Figure 6. Reaction triangle showing the advancement toward the equilibrium point for initial batch reactor concentrations of: pure A, pure B, and pure C.D. Concentration Profiles with Induced ErrorIn laboratory environments it would be understood that occasionally noise infiltrates the data obtained from the experimental measurements. Figure 7, Figure 8, and Figure 9 show the relationship of the original concentration profiles to estimated profiles with ?t=1 minute/sample and an added 5%, 10%, and 20% artificial experimental error added into the measurements respectively. It is noted that the 5% error causes the estimated reaction rates to be less accurate, and therefore the concentration profiles obtained using them are correspondingly less accurate. The estimated reaction rates for 10% error are largely more affected and made far less accurate than with 5% error. Finally 20% error makes the reaction rates so inaccurate that the concentration profiles obtained from them appear to be a completely different reacting system altogether. Note the actual reaction rates obtained in table 1 of this section for 5%, 10%, and 20% error.-635-647Figure 8. This figure shows the synthetic experimental data for concentration species with no error, then with induced artificial error of 10%, then the concentration profiles obtained with the estimated reaction rates of the error data.0Figure 8. This figure shows the synthetic experimental data for concentration species with no error, then with induced artificial error of 10%, then the concentration profiles obtained with the estimated reaction rates of the error data.31115-1905Figure 9. This figure shows the synthetic experimental data for concentration species with no error, then with induced artificial error of 20%, then the concentration profiles obtained with the estimated reaction rates of the error data.0Figure 9. This figure shows the synthetic experimental data for concentration species with no error, then with induced artificial error of 20%, then the concentration profiles obtained with the estimated reaction rates of the error data.D. Reaction Triangles with Induced ErrorSimilar to part C. of the appendix, the reaction rates data obtained in the presence of error was used to generate concentration profiles, and were then plotted in a reaction triangle. The reaction triangles for 5%, 10%, and 20% error is noted in Figure 10, Figure 11, and Figure 12 respectively. As before, it can be noted the dramatic increase in inaccuracy for 5%, 10%, and 20% error in the respective reaction triangles; from initial concentrations evolution to the equilibrium concentration itself.345064792Figure 10. Reaction triangle showing the advancement toward the equilibrium point for initial batch reactor concentrations of: pure A, pure B, and pure C as well as pure A Error, pure B Error, and pure C Error all with 5% measurement errors.0Figure 10. Reaction triangle showing the advancement toward the equilibrium point for initial batch reactor concentrations of: pure A, pure B, and pure C as well as pure A Error, pure B Error, and pure C Error all with 5% measurement errors.11214312221Figure 11. Reaction triangle showing the advancement toward the equilibrium point for initial batch reactor concentrations of: pure A, pure B, and pure C as well as pure A Error, pure B Error, and pure C Error all with 10% measurement errors.0Figure 11. Reaction triangle showing the advancement toward the equilibrium point for initial batch reactor concentrations of: pure A, pure B, and pure C as well as pure A Error, pure B Error, and pure C Error all with 10% measurement errors.492192230755Figure 12. Reaction triangle showing the advancement toward the equilibrium point for initial batch reactor concentrations of: pure A, pure B, and pure C as well as pure A Error, pure B Error, and pure C Error all with 20% measurement errors.Figure 12. Reaction triangle showing the advancement toward the equilibrium point for initial batch reactor concentrations of: pure A, pure B, and pure C as well as pure A Error, pure B Error, and pure C Error all with 20% measurement errors.E. Tables Obtained from MATLABThe tables obtained from MATLAB during experimental runs are qualified by the values obtained during execution of the program. Table 1 displays the target kinetic parameters with the non-error estimates, and then 5%, 10%, or 20% error in the measurements. The percentage of error in each error parameter is calculated by its disparity form the desired value, and the trend of the table shows that with increasing measurement error there is a proportional increase of error in the estimated kinetic parameters. The CPU or amount of time it took MATLAB to process the script also increases with the amount of measurement error. Table 2 shows the consistency of the kinetic parameters estimated with various experimental runs at specified initial concentrations. It also shows the decrease in accuracy of the estimated kinetic rate with the proportional decrease of samples taken during the reaction. The CPU varies with initial concentration and selected sampling rate.20685482287Table 1. This table shows the obtained reaction rate estimations, given the conditions of ?t=1min/sample and initial species concentrations CA0 =1, CB0=0, CC0 =0.00Table 1. This table shows the obtained reaction rate estimations, given the conditions of ?t=1min/sample and initial species concentrations CA0 =1, CB0=0, CC0 =0.207010132080kiOriginalEstimated5% ErrorPercent InaccuratekAB0.30.2920.27974.199%kac0.160.15640.15650.076%kBC0.10.09940.011788.26%kba0.10.09750.143146.83%kCA0.10.09750.012786.99%kcb0.150.14960.093337.60% kiOriginalEstimated10% ErrorPercent InaccuratekAB0.30.2920.28392.750%kac0.160.15640.15461.166%kBC0.10.0994-0.0043104.2%kba0.10.09750.067530.77%kCA0.10.09750.124727.92%kcb0.150.1496-0.0516134.4%kiOriginalEstimated20% ErrorPercent InaccuratekAB0.30.2920.249514.53%kac0.160.15640.17199.925%kBC0.10.0994-0.1619262.7%kba0.10.09750.135939.39%kCA0.10.09750.013785.93%kcb0.150.1496-0.1533202.4%00kiOriginalEstimated5% ErrorPercent InaccuratekAB0.30.2920.27974.199%kac0.160.15640.15650.076%kBC0.10.09940.011788.26%kba0.10.09750.143146.83%kCA0.10.09750.012786.99%kcb0.150.14960.093337.60% kiOriginalEstimated10% ErrorPercent InaccuratekAB0.30.2920.28392.750%kac0.160.15640.15461.166%kBC0.10.0994-0.0043104.2%kba0.10.09750.067530.77%kCA0.10.09750.124727.92%kcb0.150.1496-0.0516134.4%kiOriginalEstimated20% ErrorPercent InaccuratekAB0.30.2920.249514.53%kac0.160.15640.17199.925%kBC0.10.0994-0.1619262.7%kba0.10.09750.135939.39%kCA0.10.09750.013785.93%kcb0.150.1496-0.1533202.4%18978135105Table 2. Listed below are initial concentrations of the three species for six different experimental runs. The estimated kinetic parameters for each experimental run were evaluated for accuracy and computation time in MATLAB for sampling rates of ?t = 1, 2, and 4 minutes/sample.00Table 2. Listed below are initial concentrations of the three species for six different experimental runs. The estimated kinetic parameters for each experimental run were evaluated for accuracy and computation time in MATLAB for sampling rates of ?t = 1, 2, and 4 minutes/sample.18923081915CA0CB0CC0?tkABkackBCkbakCAkcbCPU(sec)10010.29200.1560.0990.0970.0970.1492.31001010.29200.1560.0990.0970.0970.1492.25200110.29200.1560.0990.0970.0970.1492.2280.50.5010.29200.1560.0990.0970.0970.1492.16700.50.510.29200.1560.0990.0970.0970.1492.330CA0CB0CC0?tkABkackBCkbakCAkcbCPU10020.27060.14670.09760.09070.09070.14802.289801020.27060.14670.09760.09070.09070.14802.118400120.27060.14670.09760.09070.09070.14802.29430.50.5020.27060.14670.09760.09070.09070.14802.220400.50.520.27060.14670.09760.09070.09070.14802.4353CA0CB0CC0?tkABkackBCkbakCAkcbCPU10040.21220.11940.09040.07210.07210.13972.218101040.21220.11940.09040.07210.07210.13972.151100140.21220.11940.09040.07210.07210.13972.16320.50.5040.21220.11940.09040.07210.07210.13972.253400.50.540.21220.11940.09040.07210.07210.13972.319100CA0CB0CC0?tkABkackBCkbakCAkcbCPU(sec)10010.29200.1560.0990.0970.0970.1492.31001010.29200.1560.0990.0970.0970.1492.25200110.29200.1560.0990.0970.0970.1492.2280.50.5010.29200.1560.0990.0970.0970.1492.16700.50.510.29200.1560.0990.0970.0970.1492.330CA0CB0CC0?tkABkackBCkbakCAkcbCPU10020.27060.14670.09760.09070.09070.14802.289801020.27060.14670.09760.09070.09070.14802.118400120.27060.14670.09760.09070.09070.14802.29430.50.5020.27060.14670.09760.09070.09070.14802.220400.50.520.27060.14670.09760.09070.09070.14802.4353CA0CB0CC0?tkABkackBCkbakCAkcbCPU10040.21220.11940.09040.07210.07210.13972.218101040.21220.11940.09040.07210.07210.13972.151100140.21220.11940.09040.07210.07210.13972.16320.50.5040.21220.11940.09040.07210.07210.13972.253400.50.540.21220.11940.09040.07210.07210.13972.3191MATLAB Codes Parameter_Estimation ticclear all; clc;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Target Rxn RatesK_AB=.3; K_BC=.1; K_CA=.1;K_ba=.1; K_cb=.15; K_ac=.16;%K_AB=4; K_BC=3; K_CA=4;%K_ba=6.66667; K_cb=4.5; K_ac=1.6;%Kinetic Rxn BehaviourA1 =[ -(K_AB+K_ac) K_ba K_CA ];B1 =[ K_AB -(K_ba+K_BC) K_cb ];C1 =[ K_ac K_BC -(K_CA+K_cb) ];%Generate Matrix Based on RxnsM=[A1; B1; C1];%initial values%{A0=[1 0 0 .5 .5 .1788];B0=[0 1 0 .5 0 .5054];C0=[0 0 1 0 0 .3158];%} Y0=[0.1786; 0.5051; 0.3163]; % Generate true species concentration% profile behaviourdt=.5; tf=15;[ ~,~,~,C_A,C_B,C_C,~,~,D3 ]=...ConcProf( M,Y0,dt,tf ); figure;t=0:dt:tf; plot(t,C_A,'b-*',...t,C_B,'r-*',t,C_C,'m-*');hold on;% Define Sampling Ratedt=4; tf=15;[ v1,v2,v3,C_A,C_B,C_C,D1,D2,D3 ]=...ConcProf( M,Y0,dt,tf ); %Lagrange Interpolation of 4 data pointst=0:dt:tf;[ tnew,L_A,L_B,L_C ]=... LagrangeInt_4pt( 1,tf,t,C_A,C_B,C_C );plot(tnew,L_A,'c-*',tnew,L_B,'g-*',tnew,L_C,'k-*','LineWidth',2);set(gca,'FontSize',12)title('Lagrange Interpolation of Four Data Points','FontSize',12);xlabel('Time(minutes)','FontSize',12); ylabel('Concentration','FontSize',12); hleg123 = legend('C_A(t)','C_B(t)','C_C(t)',... 'C_A(t)Lagrange','C_B(t)Lagrange','C_C(t)Langrange');hold off; figure;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Parameter Estimation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[k_AB,k_BC,k_CA,k_ba,k_ac,k_cb]=...EstRxnRates(dt,C_A,C_B,C_C);rxn_ratesGE=[k_AB,k_ac,k_BC,k_ba,k_CA,k_cb];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Estimate Concentration Behaviour %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[ C_Ae,C_Be,C_Ce ] = ConcProfEst( dt,...tf,Y0,k_AB,k_BC,k_CA,k_ba,k_ac,k_cb );Y1=Y0'; r='C_0=[ '; z=num2str(Y1);w=num2str(dt); q='] , Sampling Rate: ';L=' Minutes/Sample'; wq=[r z q w L];t=0:dt:tf; plot(t,C_Ae,'g-*',t,C_Be,...'c-*',t,C_Ce,'k-*'); title(wq);hleg1 = legend('C_A(t)','C_B(t)','C_C(t)',...'C_A(t)estimate','C_B(t)estimate',...'C_C(t)estimate');xlabel('Time(minutes)'); ylabel('Concentration'); hold off;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Parameter Estimation with Error %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%error=.05; errt=error*100;[ C_A_err,C_B_err,C_C_err ]=... Artificial_error( error,C_A,C_B,C_C );[k_AB_err,k_BC_err,k_CA_err,...k_ba_err,k_ac_err,k_cb_err]=...EstRxnRatesErr( dt,C_A_err,...C_B_err, C_C_err );[ C_AeERR,C_BeERR,C_CeERR ] =... ConcProfEst(dt,tf,Y0,k_AB_err,...k_BC_err,k_CA_err,k_ba_err,...k_ac_err,k_cb_err);figure; t=0:dt:tf; qq=num2str(errt);qqq='Effect with error of: '; q4='%';plot(t,C_A,'b-*'); hold on; plot(t,C_A_err,'r-*'); plot(t,C_AeERR,'m-*'); hleg33=legend('No Error','Error',... 'Based on Estimated Parameters from Error Data');plot(t,C_B,'b-*'); plot(t,C_B_err,'r-*');plot(t,C_BeERR,'m-*');plot(t,C_C,'b-*'); plot(t,C_C_err,'r-*'); plot(t,C_CeERR,'m-*');zz=[qqq qq q4]; title(zz); xlabel('Time(minutes)'); ylabel('Concentration');%%%%%%%%%%%%%%%%%%%%%Plot Rxn Triangle %%%%%%%%%%%%%%%%%%%%%z=length(t); figure;x=[1;0]; y=[0;1]; plot(x,y,'k'); hold on; axis([0 1 0 1]);[ f12,f22 ] = PlotRxnTri( Y0,M,dt,tf,z );plot(f12,f22,'b-*'); hold on;%Y0=[0;1;0];[ f12,f22 ] = PlotRxnTri( Y0,M,dt,tf,z );plot(f12,f22,'r-*'); hold on;%Y0=[0;0;1];[ f12,f22 ] = PlotRxnTri( Y0,M,dt,tf,z );plot(f12,f22,'m-*');hold on; Y0=[1;0;0];[ f11,f21 ] = PlotRxnTriERR( Y0,...C_A,C_B,C_C,error,z,tf,dt );plot(f11,f21,'g-*'); hold on;%Y0=[0;1;0];[ f11,f21 ] = PlotRxnTriERR( Y0,...C_A,C_B,C_C,error,z,tf,dt );plot(f11,f21,'c-*'); hold on;%Y0=[0;0;1];[ f11,f21 ] = PlotRxnTriERR( Y0,...C_A,C_B,C_C,error,z,tf,dt );plot(f11,f21,'k-*'); title('Reaction triangle Incorporating Error'); xlabel('C_A(t)'); ylabel('C_B(t)');hleg3 = legend('','Pure A','Pure B','Pure C',... 'Pure A Error','Pure B Error','Pure C Error'); hold off;% Error percentage of kinetic rates estimated with % presence of artificial errorrxn_ratesGEer=[k_AB_err,k_ac_err,k_BC_err,...k_ba_err,k_CA_err,k_cb_err];rxn_ratesGE=[k_AB,k_ac,k_BC,k_ba,k_CA,k_cb];errork=((rxn_ratesGE-rxn_ratesGEer)./rxn_ratesGE)*100;tocFunction ConcProffunction[v1,v2,v3,C_A,C_B,C_C,D1,D2,D3]=ConcProf(M,Y0,delta,T)%Find concentration equations[V,D]=eig(M);C=V\Y0; C_A=[]; C_B=[]; C_C=[]; v1=C(1)*V(:,1); v2=C(2)*V(:,2); v3=C(3)*V(:,3); D1=D(1,1); D2=D(2,2); D3=D(3,3); for t=0:delta:T y1=v1*exp(D(1,1)*t)+v2*exp(D(2,2)*t)+v3*exp(D(3,3)*t); C_A=[C_A y1(1,:)]; C_B=[C_B y1(2,:)]; C_C=[C_C y1(3,:)]; endendFunction EstRxnRatesfunction [ k_AB,k_BC,k_CA,k_ba,k_ac,k_cb ] = EstRxnRates( dt,C_A,C_B,C_C ) %UNTITLED3 Summary of this function goes here% Detailed explanation goes hereA=[]; b=[];for i=1:(length(C_A)-1) C_At1=C_A(i); C_Bt1=C_B(i); C_Ct1=C_C(i); C_At2=C_A(i+1); C_Bt2=C_B(i+1); C_Ct2=C_C(i+1); a1=[-C_At2-C_At1 -C_At2-C_At1 0 C_Bt2+C_Bt1 C_Ct2+C_Ct1 0]; a2=[C_At2+C_At1 0 -C_Bt2-C_Bt1 -C_Bt2-C_Bt1 0 C_Ct2+C_Ct1]; a3=[0 C_At2+C_At1 C_Bt2+C_Bt1 0 -C_Ct2-C_Ct1 -C_Ct2-C_Ct1]; A0=[a1; a2; a3]; A=[A; A0]; b0=[C_At2-C_At1; C_Bt2-C_Bt1; C_Ct2-C_Ct1;]; b=[b; b0]; endA=(dt/2)*A;%Manipulate to Pseudoinverse FormatAp=(A')*A; bp=(A')*b;%Gauss Eliminationrxn_ratesGE=Ap\bp; k_AB=rxn_ratesGE(1); k_BC=rxn_ratesGE(3); k_CA=rxn_ratesGE(5); k_ba=rxn_ratesGE(4); k_cb=rxn_ratesGE(6); k_ac=rxn_ratesGE(2); endFunction ConcProfEstfunction [ C_Ae,C_Be,C_Ce ] = ConcProfEst( dt,tf,Y0,k_AB,k_BC,k_CA,k_ba,k_ac,k_cb )%Kinetic Rxn BehaviourA1 =[ -(k_AB+k_ac) k_ba k_CA ];B1 =[ k_AB -(k_ba+k_BC) k_cb ];C1 =[ k_ac k_BC -(k_CA+k_cb) ];%Generate Matrix Based on RxnsM=[A1; B1; C1];%Find concentration equations[V,D]=eig(M);C=V\Y0; C_Ae=[]; C_Be=[]; C_Ce=[]; v1=C(1)*V(:,1); v2=C(2)*V(:,2); v3=C(3)*V(:,3); %Plot a Statespace******************** for t=0:dt:tf y1=v1*exp(D(1,1)*t)+v2*exp(D(2,2)*t)+v3*exp(D(3,3)*t); C_Ae=[C_Ae y1(1,:)]; C_Be=[C_Be y1(2,:)]; C_Ce=[C_Ce y1(3,:)]; endendFunction Artificial_errorfunction [C_A_err,C_B_err,C_C_err]=Artificial_error(error,C_A,C_B,C_C)a=error; b=-error;f=length(C_A);A_err=C_A.*((b-a).*rand(f,1) + a)';B_err=C_B.*((b-a).*rand(f,1) + a)';C_err=C_C.*((b-a).*rand(f,1) + a)';C_A_err= C_A+A_err;C_B_err= C_B+B_err;C_C_err= C_C+C_err;endPlotRxnTrifunction [ f12,f22 ] = PlotRxnTri( Y0,M,dt,tf,z )[ ~,~,~,C_A,C_B,C_C,~,~,~ ]=ConcProf( M,Y0,dt,tf ); [ k_AB,k_BC,k_CA,k_ba,k_ac,k_cb ] = EstRxnRates( dt,C_A,C_B,C_C );[ C_Ae,C_Be,~ ] = ConcProfEst( dt,tf,Y0,k_AB,k_BC,k_CA,k_ba,k_ac,k_cb );f12=C_Ae(1:1:z); f22=C_Be(1:1:z);endPlotRxnTriERRfunction [ f11,f21 ] = PlotRxnTriERR( Y0,C_A,C_B,C_C,error,z,tf,dt ) [ C_A_err,C_B_err,C_C_err ] = Artificial_error( error,C_A,C_B,C_C );[ k_AB_err,k_BC_err,k_CA_err,k_ba_err,k_ac_err,k_cb_err ] = EstRxnRatesErr( dt,C_A_err, C_B_err, C_C_err );[ C_AeERR,C_BeERR,~ ] = ConcProfEst( dt,tf,Y0,k_AB_err,k_BC_err,k_CA_err,k_ba_err,k_ac_err,k_cb_err );f11=C_AeERR(1:1:z); f21=C_BeERR(1:1:z); endFunction LagrangeInt_4ptfunction [ tnew,L_A,L_B,L_C ] = LagrangeInt_4pt( step,tf,time,C_A,C_B,C_C ) L_A=[]; L_B=[]; L_C=[];for k=1:step:tf L1=C_A(1)*(((k-time(2))/(time(1)-time(2))))*... (((k-time(3))/(time(1)-time(3))))*... (((k-time(4))/(time(1)-time(4)))); L2=C_A(2)*(((k-time(1))/(time(2)-time(1))))*... (((k-time(3))/(time(2)-time(3))))*... (((k-time(4))/(time(2)-time(4)))); L3=C_A(3)*(((k-time(1))/(time(3)-time(1))))*... (((k-time(2))/(time(3)-time(2))))*... (((k-time(4))/(time(3)-time(4)))); L4=C_A(4)*(((k-time(1))/(time(4)-time(1))))*... (((k-time(2))/(time(4)-time(2))))*... (((k-time(3))/(time(4)-time(3)))); L_A(k)=L1+L2+L3+L4;end for k=1:step:tf L1=C_B(1)*(((k-time(2))/(time(1)-time(2))))*... (((k-time(3))/(time(1)-time(3))))*... (((k-time(4))/(time(1)-time(4)))); L2=C_B(2)*(((k-time(1))/(time(2)-time(1))))*... (((k-time(3))/(time(2)-time(3))))*... (((k-time(4))/(time(2)-time(4)))); L3=C_B(3)*(((k-time(1))/(time(3)-time(1))))*... (((k-time(2))/(time(3)-time(2))))*... (((k-time(4))/(time(3)-time(4)))); L4=C_B(4)*(((k-time(1))/(time(4)-time(1))))*... (((k-time(2))/(time(4)-time(2))))*... (((k-time(3))/(time(4)-time(3)))); L_B(k)=L1+L2+L3+L4;end for k=1:step:tf L1=C_C(1)*(((k-time(2))/(time(1)-time(2))))*... (((k-time(3))/(time(1)-time(3))))*... (((k-time(4))/(time(1)-time(4)))); L2=C_C(2)*(((k-time(1))/(time(2)-time(1))))*... (((k-time(3))/(time(2)-time(3))))*... (((k-time(4))/(time(2)-time(4)))); L3=C_C(3)*(((k-time(1))/(time(3)-time(1))))*... (((k-time(2))/(time(3)-time(2))))*... (((k-time(4))/(time(3)-time(4)))); L4=C_C(4)*(((k-time(1))/(time(4)-time(1))))*... (((k-time(2))/(time(4)-time(2))))*... (((k-time(3))/(time(4)-time(3)))); L_C(k)=L1+L2+L3+L4;endtnew=0:step:tf-1;end ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download