Turbines: Curve Fitting and Optimization applied to HydroElectric Power BRIEF ABSTRACT This realistic problem deals with controlling the flow rate through turbines to optimize power. It is meant to be done with Lagrange multipliers. GENERAL INFORMATION FileName: Turbines Full title: Curve Fitting and Optimization applied to HydroElectric Power Last Update: 6/4/96 Developer: Aaron D. Klebanoff, Department of Mathematics, Rose-Hulman Institute of Technology, Terre Haute IN 47803 USA. Contact: 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 Background The Great Northern Paper Company operates several hydroelectric generating stations, including one on the West Branch of the Penobscot River in central Maine. Water enters the station through a pipe from a dam and can be directed toward any or all of the three turbines. The first two turbines can each accept a volume of up to 1110 cubic feet per second, and the third can accept up to 1225 cubic feet per second. The rate of flow through the ith turbine should be at least 250 cubic feet per second (otherwise it would be shut down.) The efficiency of each turbine in converting water power to electric power is a function of the rate of flow, and the efficiency is greatest about midway between the minimum and maximum rates. The level of efficiency is different for each turbine, and the rate of flow through each is constrained by the total rate of flow through the station. As an engineering consultant to the company, you are to devise a plan for distributing the total flow to the three turbines in a way that will produce the maximum amount of power at any level of flow. From Bernoulli's equation, you know that the power for each turbine is proportional to its efficiency. The efficiency of a turbine is a measure of the rate at which it can convert water power into electric power. The constant of proportionality depends mainly on the change in elevation between the headwater and the tailwater (170 feet) with an additional small term to account for mechanical energy lost due to friction. We can represent the power in the ith turbine by Wi = (170 - 1.6 10-6 QT2) hi(Qi), where Wi = power generate by turbine i (units: kilowatts/second) hi(Qi) = turbine i efficiency (units: 1/second) Qi = rate of flow through turbine i (units: cubic feet/second) QT = constant total rate of flow (units: cubic feet/second) In order to determine the efficiency of the turbines, you tested each turbine for its efficiency at various flow rates to obtain the following data. Data In[36]:= klist = Table[i, {i, 1, 11}]; Q1list = Table[250 + (1110 - 250)/10 i, {i, 0, 10}]; h1list = {9.2641, 17.693, 25.6236, 32.9053, 39.5652, 45.7038, 51.1752, 56.0342, 60.3237, 63.9635, 67.0435}; table1 = MatrixForm[ Transpose[ {Prepend[klist, "k"], Prepend[Q1list, "Q1"], Prepend[h1list, "h1"]} ] ] Out[39]//MatrixForm= k Q1 h1 1 250 9.2641 2 336 17.693 3 422 25.6236 4 508 32.9053 5 594 39.5652 6 680 45.7038 7 766 51.1752 8 852 56.0342 9 938 60.3237 10 1024 63.9635 11 1110 67.0435 In[40]:= Q2list = Q1list; h2list = {6.5659, 15.8846, 24.494, 32.405, 39.6547, 46.1733, 51.9785, 57.142, 61.5804, 65.3569, 68.4499}; table2 = MatrixForm[ Transpose[ {Prepend[klist, "k"], Prepend[Q2list, "Q2"], Prepend[h2list, "h2"]} ] ] Out[42]//MatrixForm= k Q2 h2 1 250 6.5659 2 336 15.8846 3 422 24.494 4 508 32.405 5 594 39.6547 6 680 46.1733 7 766 51.9785 8 852 57.142 9 938 61.5804 10 1024 65.3569 11 1110 68.4499 In[43]:= Q3list = Table[250 + (1225 - 250)/10 i, {i, 0, 10}]//N; h3list = {5.0419, 16.309, 26.7736, 36.5571, 45.6094, 53.8269, 61.4214, 68.2979, 74.3379, 79.7697, 84.4032}; table3 = MatrixForm[ Transpose[ {Prepend[klist, "k"], Prepend[Q3list, "Q3"], Prepend[h3list, "h3"]} ] ] Out[45]//MatrixForm= k Q3 h3 1 250. 5.0419 2 347.5 16.309 3 445. 26.7736 4 542.5 36.5571 5 640. 45.6094 6 737.5 53.8269 7 835. 61.4214 8 932.5 68.2979 9 1030. 74.3379 10 1127.5 79.7697 11 1225. 84.4032 Tasks. Plot the data of efficiency as a function of the flow rate for each of the three turbines. Recognizing that the data for the efficiency of each turbine appears to be a parabolic function of the flow rate, find ai, bi, and ci by the method of least squares so that hi(Qi) = ai + bi Qi + ci Qi2 best fits the data in the least squares sense for each turbine. In terms of QT, find the values of Q1, Q2, and Q3 that will maximize W1 + W2 + W3 subject to the constraint Q1 + Q2 + Q3 = Qi. Hint: If you plot Qi as a function of QT, it will appear linear even though your formula for it may appear quite messy. You should make a linear approximation and/or make a table of optimum flow rates as a function of QT for quick reference. Give an interval of values for QT for which your solution is valid. (You should note that for smaller values of QT, one or more of the turbines must be shut down, or for larger values of QT, one or more of the turbines must be left running indefinitely. To devise a plan to account for these extreme cases is a problem too involved to be concerned with at this time. Stick to the task of devising a plan for QT in an interval so that none of the turbines are always running and none of the turbines are shut down.) Write a report to the engineer in charge of the project detailing your recommendations. Leave details for an appendix. KEYWORDS Lagrange Multipliers, Least Squares, Curve Fitting, Optimization, Power. TEACHER NOTES ISSUES RELATED TO THE PROBLEM Prerequisites This problem works best if the students have already been taught how to use Lagrange multipliers. Students must be familiar with the method of least squares. Students will find the curve fitting part of this project much easier if in addition to the usual linear regression, you also fit data to a quadratic or cubic polynomial. Although the problem was intended to have a curve fitting component, curve fitting software could be used to save time. Time allotment - time management Allow 1 - 2 in class days + 1 - 2 weeks out of class for this project. On the first day, allow the students to become familiar with the problem. It is this day that students will listen attentively to a demonstration on how to use least squares for fitting curves other than lines. If you have time to devote 2 class periods for this project, the 2nd one should be after the students have had time to work on the project on their own so that they can check with you on their progress as well as see how the other teams are doing. Expectations Students like the realism of the problem. Students will take pride in their final reports. Some students may require additional help to understand the problem. In particular, some will have difficulty with the notation, some will not recognize that the tasks are designed to walk the student through the solution. Finally, some students will be confused by the fourth task in which they are to determine the range of river flow rates for which their solution is valid. All of these issues are relatively easy to clear up when explained in "real-life" terms. In summary, it takes time for students to become familiar enough with the problem so that they understand what their "tasks" are for. Future payoffs Students learn that the method of least squares is not just for doing linear fits. Students see how several tools can be used together to solve a real-life problem. Students gain experience writing business reports. Extensions Bright students should attempt to solve the problem for any river flow rate. It's important that they consider examples to show why they cannot simply extend their solutions to the original problem. References and Sources Calculus, 2nd ed. by Richard A. Hunt, 1994. Harper Collins, New York, p. 881. POSSIBLE SOLUTION Data From Project In[46]:= x1 = Table[250 + (1110 - 250)/10 i, {i, 0, 10}] Out[46]= {250, 336, 422, 508, 594, 680, 766, 852, 938, 1024, 1110} In[47]:= x2 = Table[250 + (1110 - 250)/10 i, {i, 0, 10}] Out[47]= {250, 336, 422, 508, 594, 680, 766, 852, 938, 1024, 1110} In[48]:= x3 = Table[250 + (1225 - 250)/10 i, {i, 0, 10}]//N Out[48]= {250., 347.5, 445., 542.5, 640., 737.5, 835., 932.5, 1030., 1127.5, 1225.} In[49]:= y1 = {9.2641, 17.693, 25.6236, 32.9053, 39.5652, 45.7038, 51.1752, 56.0342, 60.3237, 63.9635, 67.0435}; In[50]:= y2 = {6.5659, 15.8846, 24.494, 32.405, 39.6547, 46.1733, 51.9785, 57.142, 61.5804, 65.3569, 68.4499}; In[51]:= y3 = {5.0419, 16.309, 26.7736, 36.5571, 45.6094, 53.8269, 61.4214, 68.2979, 74.3379, 79.7697, 84.4032}; Least Square Fits of the Data For h1(q1) In[52]:= n = 11; x = x1; y = y1; In[55]:= xsum = Sum[x[[i]], {i, 1, n}]; ysum = Sum[y[[i]], {i, 1, n}]; x2sum = Sum[x[[i]]^2, {i, 1, n}]; x3sum = Sum[x[[i]]^3, {i, 1, n}]; x4sum = Sum[x[[i]]^4, {i, 1, n}]; xysum = Sum[x[[i]] y[[i]], {i, 1, n}]; x2ysum = Sum[x[[i]]^2 y[[i]], {i, 1, n}]; In[62]:= rel1 = c x4sum + b x3sum + a x2sum - x2ysum; rel2 = c x3sum + b x2sum + a xsum - xysum; rel3 = c x2sum + b xsum + a n - ysum; In[65]:= eqn1 = rel1 == 0; eqn2 = rel2 == 0; eqn3 = rel3 == 0; In[68]:= coefs = Flatten[Solve[{eqn1, eqn2, eqn3}, {a, b, c}]] Out[68]= {a -> -18.8874, b -> 0.122675, c -> -0.0000407725} In[69]:= a1 = a /. coefs; b1 = b /. coefs; c1 = c /. coefs; In[72]:= nh1 = a1 + b1 q1 + c1 q1^2; In[73]:= points1 = Table[{x1[[i]], y1[[i]]}, {i, 1, n}]; In[74]:= pts1 = ListPlot[points1] Out[74]= -Graphics- In[75]:= curve1 = Plot[nh1, {q1, x1[[1]], x1[[n]]}] Out[75]= -Graphics- In[76]:= Show[pts1, curve1] Out[76]= -Graphics- For h2(q2) In[77]:= n = 11; x = x2; y = y2; In[80]:= xsum = Sum[x[[i]], {i, 1, n}]; ysum = Sum[y[[i]], {i, 1, n}]; x2sum = Sum[x[[i]]^2, {i, 1, n}]; x3sum = Sum[x[[i]]^3, {i, 1, n}]; x4sum = Sum[x[[i]]^4, {i, 1, n}]; xysum = Sum[x[[i]] y[[i]], {i, 1, n}]; x2ysum = Sum[x[[i]]^2 y[[i]], {i, 1, n}]; In[87]:= rel1 = c x4sum + b x3sum + a x2sum - x2ysum; rel2 = c x3sum + b x2sum + a xsum - xysum; rel3 = c x2sum + b xsum + a n - ysum; In[90]:= eqn1 = rel1 == 0; eqn2 = rel2 == 0; eqn3 = rel3 == 0; In[93]:= coefs = Flatten[Solve[{eqn1, eqn2, eqn3}, {a, b, c}]] Out[93]= {a -> -24.3935, b -> 0.135594, c -> -0.000046819} In[94]:= a2 = a /. coefs; b2 = b /. coefs; c2 = c /. coefs; In[97]:= nh2 = a2 + b2 q2 + c2 q2^2; In[98]:= points2 = Table[{x2[[i]], y2[[i]]}, {i, 1, n}]; In[99]:= pts2 = ListPlot[points2] Out[99]= -Graphics- In[100]:= curve2 = Plot[nh2, {q2, x2[[1]], x2[[n]]}] Out[100]= -Graphics- In[101]:= Show[pts2, curve2] Out[101]= -Graphics- For h3(q3) In[102]:= n = 11; x = x3; y = y3; In[105]:= xsum = Sum[x[[i]], {i, 1, n}]; ysum = Sum[y[[i]], {i, 1, n}]; x2sum = Sum[x[[i]]^2, {i, 1, n}]; x3sum = Sum[x[[i]]^3, {i, 1, n}]; x4sum = Sum[x[[i]]^4, {i, 1, n}]; xysum = Sum[x[[i]] y[[i]], {i, 1, n}]; x2ysum = Sum[x[[i]]^2 y[[i]], {i, 1, n}]; In[112]:= rel1 = c x4sum + b x3sum + a x2sum - x2ysum; rel2 = c x3sum + b x2sum + a xsum - xysum; rel3 = c x2sum + b xsum + a n - ysum; In[115]:= eqn1 = rel1 == 0; eqn2 = rel2 == 0; eqn3 = rel3 == 0; In[118]:= coefs = Flatten[Solve[{eqn1, eqn2, eqn3}, {a, b, c}]] Solve::luc: -- Message text not found -- (Solve) ( 8 -4.75448 10 + <<3>> == 0 && <<2>>) Out[118]= {a -> -27.0525, b -> 0.138083, c -> -0.0000384533} In[119]:= a3 = a /. coefs; b3 = b /. coefs; c3 = c /. coefs; In[122]:= nh3 = a3 + b3 q3 + c3 q3^2; In[123]:= points3 = Table[{x3[[i]], y3[[i]]}, {i, 1, n}]; In[124]:= pts3 = ListPlot[points3] Out[124]= -Graphics- In[125]:= curve3 = Plot[nh3, {q3, x3[[1]], x3[[n]]}] Out[125]= -Graphics- In[126]:= Show[pts3, curve3] Out[126]= -Graphics- Maximize Power with Lagrange Multipliers Maximize power function, powtot: In[127]:= dh = 170; e = 1.6 10^(-6); In[129]:= pow1 = (dh - e qt^2) nh1; pow2 = (dh - e qt^2) nh2; pow3 = (dh - e qt^2) nh3; powtot = pow1 + pow2 + pow3; Subject to the condition that q1 + q2 + q3 = qt; Write as g = 0. In[133]:= g = q1 + q2 + q3 - qt; Do this by minimizing Lagrange function, lag: In[134]:= lag = powtot - lambda g; Set derivatives equal to zero: Derivatives: In[135]:= dlq1 = D[lag, q1]; In[136]:= dlq2 = D[lag, q2]; In[137]:= dlq3 = D[lag, q3]; In[138]:= dllam = D[lag, lambda]; Equations: In[139]:= eq1 = dlq1 == 0; In[140]:= eq2 = dlq2 == 0; In[141]:= eq3 = dlq3 == 0; In[142]:= eq4 = dllam == 0; Solve For the Critical Point(s): In[143]:= cp = Flatten[Simplify[ Solve[{eq1, eq2, eq3, eq4}, {q1, q2, q3, lambda}] ]] Out[143]= {q1 -> (0.341161 (-10307.8 + qt) (-10307.8 + qt) (-10307.8 + qt) (-338.308 + qt) (10307.8 + qt) (10307.8 + qt) (10307.8 + qt)) / 24 6 (-1.19946 10 + 1.17965 10 qt + 16 2 3 3.38672 10 qt - 0.0100098 qt - 8 4 6 3.1875 10 qt + qt ), q2 -> 1448.07 - (1. -6 (-0.0136982 + 2.88509 10 qt + -10 2 -14 3 3.86773 10 qt - 8.14614 10 qt - -18 4 -22 5 3.64021 10 qt + 7.66695 10 qt + -26 6 -30 7 1.14203 10 qt - 2.40532 10 qt )) / -6 -13 2 (-9.71078 10 + 2.74187 10 qt - -21 4 -30 6 2.58058 10 qt + 8.09594 10 qt ), q3 -> 1795.47 - (1. -6 (-0.0136982 + 2.88509 10 qt + -10 2 -14 3 3.86773 10 qt - 8.14614 10 qt - -18 4 -22 5 3.64021 10 qt + 7.66695 10 qt + -26 6 -30 7 1.14203 10 qt - 2.40532 10 qt )) / -6 -13 2 (-7.97565 10 + 2.25195 10 qt - -21 4 -30 6 2.11948 10 qt + 6.64935 10 qt ), lambda -> -6 (-1. (-0.0136982 + 2.88509 10 qt + -10 2 -14 3 3.86773 10 qt - 8.14614 10 qt - -18 4 -22 5 3.64021 10 qt + 7.66695 10 qt + -26 6 -30 7 1.14203 10 qt - 2.40532 10 qt )) / -11 2 (0.000610034 - 1.1483 10 qt + -20 4 5.40376 10 qt )} In[144]:= q1max = Factor[q1 /. cp] Out[144]= (0.341161 (-10307.8 + qt) (-10307.8 + qt) (-10307.8 + qt) (-338.308 + qt) (10307.8 + qt) (10307.8 + qt) (10307.8 + qt)) / ((-10307.8 + qt) (-10307.8 + qt) (-10307.8 + qt) (10307.8 + qt) (10307.8 + qt) (10307.8 + qt)) In[145]:= q1maxap = Normal[Series[q1max, {qt, 0, 1}]] Out[145]= -115.417 + 0.341161 qt In[146]:= q2max = Factor[q2 /. cp] Out[146]= (0.297102 (-10307.8 + qt) (-10307.8 + qt) (-10307.8 + qt) (126.057 + qt) (10307.8 + qt) (10307.8 + qt) (10307.8 + qt)) / ((-10307.8 + qt) (-10307.8 + qt) (-10307.8 + qt) (10307.8 + qt) (10307.8 + qt) (10307.8 + qt)) In[147]:= q2maxap = Normal[Series[q2max, {qt, 0, 1}]] Out[147]= 37.4518 + 0.297102 qt In[148]:= q3max = Factor[q3 /. cp] Out[148]= (0.361737 (-10307.8 + qt) (-10307.8 + qt) (-10307.8 + qt) (215.531 + qt) (10307.8 + qt) (10307.8 + qt) (10307.8 + qt)) / ((-10307.8 + qt) (-10307.8 + qt) (-10307.8 + qt) (10307.8 + qt) (10307.8 + qt) (10307.8 + qt)) In[149]:= q3maxap = Normal[Series[q3max, {qt, 0, 1}]] Out[149]= 77.9656 + 0.361737 qt Graphs of the optimum flow rates qi as functions of the river flow rate qt. In[150]:= Plot[q1max, {qt, 750, 3445}, PlotRange -> {250, 1110}] Out[150]= -Graphics- In[151]:= Plot[q2max, {qt, 750, 3445}, PlotRange -> {250, 1110}] Out[151]= -Graphics- In[152]:= Plot[q3max, {qt, 750, 3445}, PlotRange -> {250, 1225}] Out[152]= -Graphics- The following is just a check that all water is being used... In[153]:= Plot[q1max + q2max + q3max, {qt, 750, 3445}] Out[153]= -Graphics- Intervals On Which All Three are at partial capacity. Minimum Input Flow Rates for q1. In[154]:= FindRoot[q1max == 250, {qt, 1050}] Out[154]= {qt -> 1071.1} In[155]:= q1maxlow = Max[qt /. %, 3 250] Out[155]= 1071.1 Maximum Input Flow Rates for q1. In[156]:= FindRoot[q1max == 1110, {qt, 1050}] Out[156]= {qt -> 3591.9} In[157]:= q1maxhigh = Min[qt /. %, 1110 + 1110 + 1225] Out[157]= 3445 Minimum Input Flow Rates for q2. In[158]:= FindRoot[q2max == 250, {qt, 1050}] Out[158]= {qt -> 715.406} In[159]:= q2maxlow = Max[qt /. %, 3 250] Out[159]= 750 Maximum Input Flow Rates for q2. In[160]:= FindRoot[q2max == 1110, {qt, 1050}] Out[160]= {qt -> 3610.04} In[161]:= q2maxhigh = Min[qt /. %, 1110 + 1110 + 1225] Out[161]= 3445 Minimum Input Flow Rates for q3. In[162]:= FindRoot[q3max == 250, {qt, 1050}] Out[162]= {qt -> 475.578} In[163]:= q3maxlow = Max[qt /. %, 3 250] Out[163]= 750 Maximum Input Flow Rates for q3. In[164]:= FindRoot[q3max == 1225, {qt, 1050}] Out[164]= {qt -> 3170.9} In[165]:= q3maxhigh = Min[qt /. %, 1110 + 1110 + 1225] Out[165]= 3170.9 Minimum and Maximum Input Flow Rates for q1, q2, and q3. In[166]:= qtmin = Max[q1maxlow, q2maxlow, q3maxlow] Out[166]= 1071.1 In[167]:= qtmax = Max[q1maxhigh, q2maxhigh, q3maxhigh] Out[167]= 3445 ISSUES IN SOLUTION There are many ways to do the linear approximation. Most students will choose two points on their curve (which looks like a line) --- hopefully at end points of the feasible interval --- to get their line. There are at least three common ways to present the solution to the head engineer: Include (1) graphical information, (2) tabular information, or (3) analytical (i.e. formulas) information. Many students will choose a combination of these methods. Allow the students to decide what they believe is the most appropriate way to present their findings. You might suggest that the head engineer may want to program the formulas into a computer system so it will be necessary to add these formulas in an appendix of the report.