CHEMOPT BRIEF ABSTRACT After modeling a chemical reaction A->B->C with a set of linear, first-order differential equations we estimate kinetic parameters through a number of methods and determine the run time for the reaction which will optimize a financial return on the reaction. GENERAL INFORMATION FileName: CHEMOPT Full title: Reaction chemistry, parameter estimation, and optimization Last Revision Date: 3 June 1996. Developer: Brian J. Winkel, Department of Mathematical Sciences, United States Military Academy, West Point NY 10996 USA. Phone: 914-938-3200. Email: ab3646@usma2.usma.edu. FAX: 914-938-2409. Contact: Brian J. Winkel, Department of Mathematical Sciences, United States Military Academy, West Point NY 10996 USA. Phone: 914-938-3200. Email: ab3646@usma2.usma.edu. FAX: 914-938-2409. Aaron D. Klebanoff, Department of Mathematics, Rose-Hulman Institute of Technology, Terre Haute IN 47803 USA. Phone: 812-877-8151. Email: Aaron.Klebanoff@Rose-Hulman.Edu. FAX: 812-877-3198. Support: The production of this material is supported by the National Science Foundation under Division of Undergraduate Education grant DUE-9352849: Development Site for Complex, Technology-Based Problems in Calculus with Applications in Science and Engineering and the Arvin Foundation of Columbus IN. STATEMENT OF PROBLEM A laboratory experiment is going on in the Projects Lab of your company. A colleague, a production chemist, comes to you for advice. Compound A is heated to 120oC in order to prodcue compounds B and C. The temperature of the pot containing all of these compounds is kept at 120oC in order to keep the reaction going. It is believed the reaction is simple first order reaction where the reaction rate of converting compound A to B is k1 and the reaction rate of converting compound B to C is k2. The following data has been gathered Time (min) A (mole) B (mole) C(mole) 0 1.00 0.00 0.000 2 0.88 0.12 0.003 6 0.69 0.29 0.030 10 0.53 0.42 0.050 20 0.28 0.56 0.16 30 0.15 0.57 0.28 50 0.043 0.46 0.50 70 0.012 0.33 0.66 90 0.22 0.78 120 0.12 0.88 150 0.06 0.94 200 0.02 0.98 The goal is to produce compound B for marketing, thus we seek a mathematical model as an aid in telling us the best time to stop the process and extract compound B for the market. But a measure of best is most profit or net return. We know that A costs $0.50 per mole and B sells for $3.50 per mole. Our production engineer friend has found that for compound C we can expect a return of $0.25 per mole from the recovered C. Find the optimal shutoff time for this process in order to optimize the net return on this process. Suggestions: (1) Without using any "sophisticated" modelling techniques offer a plausible approach and suggest a method to solve such a problem. Look at the data! (2) Form a mathematical model. Be sure to state your assumptions - mathematical and chemical. Use this mathematical model to tell your corporate friend just when to shut off the process and just how much profit they can expect on this process. Ideas for modeling consist of curve fitting, differential (rate) equations, or difference equations. Post Solution issues, modifications, and more problems: Suppose your colleague comes back to you and says that you NOW have to consider the optimnization in view of the energy costs to run the process and it is known that the energy costs to run this process are $0.005/min. Formulate a model which will find the optimal shutoff time for this process in order to optimize the net return on this process. Compare your results when you consider energy to when you do not. Finally attempt to determine the sensitivity of yield and final net return in terms of changes in the energy cost. Namely offer a plot of of net return and change in net return/unit change in price of energy (marginal net return) as a function of price of energy, say in the range [0, .040], i.e. from no costs to approximately 8 times the current costs of $0.005/min. KEYWORDS Chemical reaction, first-order kinetics, rate (differential, difference) equation, data fitting, optimization, sensitivity. TEACHER NOTES ISSUES RELATED TO THE PROBLEM Prerequisites Some exposure to some notion that the rate of reaction of a decaying compound is proportional to the amount of the compound present, perhaps (but not essential) Law of Mass Action or first order kinetics in chemistry course. Some exposure to the notion of a rate or differential equation and a system of differential equations. Use DSolve (or other CAS differential equation solver) for the system of differential equations and thus an understanding of what a "solution to a differential equation or a system of differential equations" means. Some understanding of the concept of least squares fitting of functions, nonlinear in the parameters to be determined, to data. Optimization of a function of one variable. Time allotment - time management Students can cvonstruct the system of differential equations in 10 minutes if they have been coached to think about Law of Mass Action or rate equations. The analysis of solving (using DSolve say in Mathematica) and the least squares analysis needs introduction, but if they have had some experience with curve fitting to data this material could be introduced in one hour and bascially assigned for homework over a few nights. Expectations We would expect students to formulate a model of the amounts or reactions for the chemicals involved and then estimate the parameters in the model. Future payoffs Students will be able to formulate models in which unknown parameters exist and these parameters may be found using data. Extensions Work on large system of chemicals or make costs nonlinear, i.e., cheaper in larger lots. References and Sources POSSIBLE SOLUTION(S) Idea - not using differential equation model, but direct curve fitting. We turn to this laters. We could solve just the differential equation for a(t) - exponential decay - and then use this with antidifferentiation (e.g., integrating factor) to get b(t) and then easily get c(t). We set up differential or rate equations for the chemical reactions assuming rate of reaction is proportional to amount of chemical present. We assume the reactions are first order. eqa = a'[t] == - k1 a[t] a'[t] == -(k1 a[t]) a'[t] == -(k1 a[t]) a'[t] == -(k1 a[t]) eqb = b'[t] == - k2 b[t] + k1 a[t] b'[t] == k1 a[t] - k2 b[t] b'[t] == k1 a[t] - k2 b[t] b'[t] == k1 a[t] - k2 b[t] eqc = c'[t] == k2 b[t] c'[t] == k2 b[t] c'[t] == k2 b[t] c'[t] == k2 b[t] We solve the differential equations for a, b, and c as functions of time with the rate constants k1 and k2 in each. sol = DSolve[{eqa,eqb,eqc,a[0]==1,b[0]==0, c[0]==0}, {a[t],b[t],c[t]},t] 1 k1 {{a[t] -> -------------- - -----------------, k1 t k1 k1 t k1 E (1 - --) E (1 - --) k2 k2 k2 k1 k1 b[t] -> ----------------- - -----------------, k1 t k1 k2 t k1 E (1 - --) k2 E (1 - --) k2 k2 k2 1 k1 c[t] -> 1 - -------------- + -----------------}} k1 t k1 k2 t k1 E (1 - --) E (1 - --) k2 k2 k2 We grab off each solution: asol[t_] = Simplify[a[t]/.sol[[1]]]; bsol[t_] = Simplify[b[t]/.sol[[1]]]; csol[t_] = Simplify[c[t]/.sol[[1]]]; We input and plot all the data. adata = {{0,1},{2,.88},{6,.69},{10,.53}, {20,.28},{30,.15},{50,.043},{70,.012}, {90,0},{120,0},{150,0},{200,0}}; bdata = {{0,1},{2,.12},{6,.29},{10,.42}, {20,.55},{30,.57},{50,.46},{70,.33}, {90,.22},{120,.12},{150,.06},{200,.02}}; cdata = {{0,0},{2,.003},{6,.030},{10,.050}, {20,.16},{30,.28},{50,.50},{70,.66}, {90,.78},{120,.88},{150,.94},{200,.98}}; adataPlot = ListPlot[adata,PlotStyle->PointSize[.01]] -Graphics- bdataPlot = ListPlot[bdata,PlotStyle->PointSize[.015]] -Graphics- cdataPlot = ListPlot[cdata,PlotStyle->PointSize[.02]] -Graphics- alldataPlot =Show[adataPlot,bdataPlot,cdataPlot] -Graphics- We use the data to determine the rate constants k1 and k2. We formulate a least squares function in the constants (now variables) k1 and k2 using the data from the chemicals -- first separately and then together. We formulate a least squares function in the constants (now variables) k1 and k2 using the data from the chemical a. SSa[k1_,k2_] = Sum[(asol[adata[[n,1]]] - adata[[n,2]])^2, {n,1, Length[adata]}]; solak = FindMinimum[SSa[k1,k2],{k1,.1},{k2,.2}] {0.0000535918, {k1 -> 0.063221, k2 -> 0.2}} We formulate a least squares function in the constants (now variables) k1 and k2 using the data from the chemical b. SSb[k1_,k2_] = Sum[(bsol[bdata[[n,1]]] - bdata[[n,2]])^2, {n,1, Length[bdata]}]; solbk = FindMinimum[SSb[k1,k2],{k1,.1},{k2,.2}] {1.00015, {k1 -> 0.0623575, k2 -> 0.0210941}} We formulate a least squares function in the constants (now variables) k1 and k2 using the data from the chemical c. SSc[k1_,k2_] = Sum[(csol[cdata[[n,1]]] - cdata[[n,2]])^2, {n,1, Length[cdata]}]; solck = FindMinimum[SSc[k1,k2],{k1,.1},{k2,.2}] {0.000135185, {k1 -> 0.0210623, k2 -> 0.0640486}} Now we formulate a least squares function in the constants (now variables) k1 and k2 using the data from ALL the chemicals. SS[k1_,k2_] = SSa[k1,k2] + SSb[k1,k2] + SSc[k1,k2]; solk = FindMinimum[SS[k1,k2],{k1,.1},{k2,.2}] {1.00038, {k1 -> 0.0629595, k2 -> 0.0211523}} The rate constants k1 and k2 seem to differ depending upon which of the data sets we wish to "satisfy." Thus we take the values which satisfy "all" of the data. We plug the values for the rate constants into the solved a, b, and c functions. cfit[t_] = csol[t]/.solk[[2]]; bfit[t_] = bsol[t]/.solk[[2]]; afit[t_] = asol[t]/.solk[[2]]; And plot the solutions - over the data. fPlot = Plot[{afit[t],bfit[t],cfit[t]},{t,0,200}] -Graphics- Show[alldataPlot,fPlot] -Graphics- This indicates we have a good model for the fit is almost perfect. We optimize the objective function - no energy costs. P[t_] = -.5 asol[t] + 3.5 bsol[t] + .25 csol[t]/.solk[[2]]; Plot[P[t],{t,0,200}] -Graphics- solp = FindRoot[P'[t]==0,{t,40}] {t -> 29.5001} P[t]/.solp[[1]] 1.99135 We optimize the objective function - with energy costs. Pcost[t_] = -.005 t -.5 asol[t] + 3.5 bsol[t] + .25 csol[t]/.solk[[2]]; Plot[Pcost[t],{t,0,200}] -Graphics- solPcost = FindRoot[Pcost'[t]==0,{t,40}] {t -> 27.5171} Pcost[t]/.solPcost[[1]] 1.84894 We now determine the sensitivity of the operation to energy costs. We compute the net return as a function of time and energy rate s $/min to run the reaction. net[t_,s_] = - s t -.5 asol[t] + 3.5 bsol[t] + .25 csol[t]/.solk[[2]]; We compute the change in net rate as a function of time. dnet[t_,s_] = D[net[t,s],t]; setStuff = {}; Do[setStuff = Append[setStuff, {.0001 i,net[t,.0001 i]}/.FindRoot[dnet[t,.0001 i]==0,{t,2,4}][[1]] ], {i,0,400,10}] setPlot = ListPlot[setStuff,PlotStyle->PointSize[.01], PlotRange->{{0,.05},{0,2.0}}, AxesLabel->{"Energy cost ($/min)","MaxNet"}] -Graphics- d[x_] = Fit[setStuff,{1,x,x^2},x] 2 1.98833 - 28.4241 x + 130.896 x dPlot = Plot[d[x],{x,0,.05}] -Graphics- Show[setPlot,dPlot] -Graphics- Now we compute the difference in net return as a function of the energy price in $/min - this is the marginal change in net return as a function of energy. diffsetStuff = {}; Do[diffsetStuff = Append[diffsetStuff, {setStuff[[i]][[1]], (setStuff[[i+1]][[2]] - setStuff[[i]][[2]])* 1/(setStuff[[i+1]][[1]] - setStuff[[i]][[1]])} ], {i,1, Length[setStuff]-1}] df[x_] = Fit[diffsetStuff,{1,x,x^2},x] 2 -29.1517 + 370.934 x - 2735.84 x dfPlot = Plot[df[x],{x,0,.05}] -Graphics- diffsetPlot = ListPlot[diffsetStuff,PlotStyle->PointSize[.01], PlotRange->{{0,.05},{-30,-18}}, AxesLabel->{"Energy cost ($/min)","MaxNet/($/min)"}] -Graphics- Show[diffsetPlot,dfPlot] -Graphics- There is a greater marginal loss in revenue when energy costs are low. Another approach from start - We could fit polynomials to the data fa[t_] = Fit[adata,{1,t,t^2,t^3,t^4,t^5},t] 2 0.983247 - 0.0538534 t + 0.00115793 t - 3 -8 4 -10 5 0.0000118174 t + 5.66228 10 t - 1.01836 10 t faPlot = Plot[fa[t],{t,0,200},PlotRange->{{0,200},{0,1}}] -Graphics- Show[faPlot,adataPlot] -Graphics- fb[t_] = Fit[bdata,{1,t,t^2,t^3,t^4},t] 2 0.491081 + 0.00114856 t - 0.0000633133 t + -7 3 -10 4 2.55029 10 t - 1.27333 10 t fbPlot = Plot[fb[t],{t,0,200},PlotRange->{{0,200},{0,1}}] -Graphics- Show[fbPlot,bdataPlot] -Graphics- fc[t_] = Fit[cdata,{1,t,t^2,t^3,t^4},t] 2 -0.0195834 + 0.0087579 t + 0.0000622727 t - -7 3 -9 4 8.81846 10 t + 2.38368 10 t fcPlot = Plot[fc[t],{t,0,200},PlotRange->{{0,200},{0,1}}] -Graphics- Show[fcPlot,cdataPlot] -Graphics- We compute a net return function based on our polynomial fitted curves. nPcost[t_] = - .005 t -.5 fa[t] + 3.5 fb[t] + .25 fc[t]; Plot[nPcost[t],{t,0,200}] -Graphics- We recall the fitted model to our differential equation solution, Pcost(t). Plot[Pcost[t],{t,0,200}] solPcost = FindRoot[Pcost'[t]==0,{t,40}] {t -> 27.5171} Pcost[t]/.solPcost[[1]] 1.84894 We now compute the optimal net return using our polynomial fitted functions and see that we do not do as well. solnPcost = FindRoot[nPcost'[t]==0,{t,40}] {t -> 24.4458} nPcost[t]/.solnPcost[[1]] 1.52749 We note that we did better, i.e., our net return was higher with the differential equation model. ISSUES IN SOLUTION Somehow students need to realize that they can formulate a mathematical model with parameters and use data to determine these parameters.