One dimensional optimization problem

one dimensional optimization algorithm and also one-dimensional cutting stock optimization in consecutive time periods
GregDeamons Profile Pic
GregDeamons,New Zealand,Professional
Published Date:03-08-2017
Your Website URL(Optional)
cha01102_ch07_182-204.qxd 12/18/10 1:57 PM Page 182 7 Optimization CHAPTER OBJECTIVES The primary objective of this chapter is to introduce you to how optimization can be used to determine minima and maxima of both one-dimensional and multidimensional functions. Specific objectives and topics covered are Understanding why and where optimization occurs in engineering and scientific • problem solving. Recognizing the difference between one-dimensional and multidimensional • optimization. Distinguishing between global and local optima. • Knowing how to recast a maximization problem so that it can be solved with a • minimizing algorithm. Being able to define the golden ratio and understand why it makes one- • dimensional optimization efficient. Locating the optimum of a single-variable function with the golden-section search. • Locating the optimum of a single-variable function with parabolic interpolation. • Knowing how to apply the fminbnd function to determine the minimum of a • one-dimensional function. Being able to develop MATLAB contour and surface plots to visualize two- • dimensional functions. Knowing how to apply the fminsearch function to determine the minimum of a • multidimensional function. YOU’VE GOT A PROBLEM n object like a bungee jumper can be projected upward at a specified velocity. If it is subject to linear drag, its altitude as a function of time can be computed as A     m mg mg −(c/m)t z = z + v + 1 − e − t (7.1) 0 0 c c c 182cha01102_ch07_182-204.qxd 12/20/10 7:21 AM Page 183 7.1 INTRODUCTION AND BACKGROUND 183 z, m 200 Maximum 100 elevation t, s 0 2468 10 12 –100 FIGURE 7.1 Elevation as a function of time for an object initially projected upward with an initial velocity. where z = altitude (m) above the earth’s surface (defined as z = 0), z = the initial altitude 0 (m), m = mass (kg), c = a linear drag coefficient (kg/s), v = initial velocity (m/s), and t = 0 time (s). Note that for this formulation, positive velocity is considered to be in the upward 2 direction. Given the following parameter values: g = 9.81 m/s , z = 100 m, v = 55 m/s, 0 0 m = 80 kg, and c = 15 kg/s, Eq. (7.1) can be used to calculate the jumper’s altitude. As displayed in Fig. 7.1, the jumper rises to a peak elevation of about 190 m at about t = 4 s. Suppose that you are given the job of determining the exact time of the peak elevation. The determination of such extreme values is referred to as optimization. This chapter will introduce you to how the computer is used to make such determinations. 7.1 INTRODUCTION AND BACKGROUND In the most general sense, optimization is the process of creating something that is as effective as possible. As engineers, we must continuously design devices and products that perform tasks in an efficient fashion for the least cost. Thus, engineers are always con- fronting optimization problems that attempt to balance performance and limitations. In addition, scientists have interest in optimal phenomena ranging from the peak elevation of projectiles to the minimum free energy. From a mathematical perspective, optimization deals with finding the maxima and minima of a function that depends on one or more variables. The goal is to determine the values of the variables that yield maxima or minima for the function. These can then be substituted back into the function to compute its optimal values. Although these solutions can sometimes be obtained analytically, most practical optimization problems require numerical, computer solutions. From a numerical stand- point, optimization is similar in spirit to the root-location methods we just covered in Chaps. 5 and 6. That is, both involve guessing and searching for a point on a function. The fundamental difference between the two types of problems is illustrated in Fig. 7.2. Root location involves searching for the location where the function equals zero. In contrast, optimization involves searching for the function’s extreme points.cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 184 184 OPTIMIZATION f(x) = 0 f (x) Maximum f(x) 0 f (x) = 0 Root 0 x Root Root f(x) = 0 Minimum f(x) 0 FIGURE 7.2 A function of a single variable illustrating the difference between roots and optima. As can be seen in Fig. 7.2, the optimums are the points where the curve is flat. In math-  ematical terms, this corresponds to the x value where the derivative f (x) is equal to zero.  Additionally, the second derivative, f (x), indicates whether the optimum is a minimum or   a maximum: if f (x) 0, the point is a maximum; if f (x) 0, the point is a minimum. Now, understanding the relationship between roots and optima would suggest a possi- ble strategy for finding the latter. That is, you can differentiate the function and locate the root (i.e., the zero) of the new function. In fact, some optimization methods do just this by  solving the root problem: f (x) = 0. EXAMPLE 7.1 Determining the Optimum Analytically by Root Location Problem Statement. Determine the time and magnitude of the peak elevation based on 2 Eq. (7.1). Use the following parameter values for your calculation: g = 9.81 m/s , z = 100 m, v = 55 m/s, m = 80 kg, and c = 15 kg/s. 0 0 Solution. Equation (7.1) can be differentiated to give   dz mg −(c/m)t −(c/m)t = v e − 1 − e (E7.1.1) 0 dt c Note that because v = dz/dt, this is actually the equation for the velocity. The maximum elevation occurs at the value of t that drives this equation to zero. Thus, the problem amounts to determining the root. For this case, this can be accomplished by setting the de- rivative to zero and solving Eq. (E7.1.1) analytically for   m cv 0 t = ln 1 + c mg Substituting the parameters gives   80 15(55) t = ln 1 + = 3.83166 s 15 80(9.81)cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 185 7.1 INTRODUCTION AND BACKGROUND 185 This value along with the parameters can then be substituted into Eq. (7.1) to compute the maximum elevation as     80 80(9.81) 80(9.81) −(15/80)3.83166 z = 100+ 50 + 1 − e − (3.83166) = 192.8609 m 15 15 15 We can verify that the result is a maximum by differentiating Eq. (E7.1.1) to obtain the second derivative 2 d z c m −(c/m)t −(c/m)t =− v e − ge =−9.81 0 2 2 dt m s The fact that the second derivative is negative tells us that we have a maximum. Further, the result makes physical sense since the acceleration should be solely equal to the force of gravity at the maximum when the vertical velocity (and hence drag) is zero. Although an analytical solution was possible for this case, we could have obtained the same result using the root-location methods described in Chaps. 5 and 6. This will be left as a homework exercise. Although it is certainly possible to approach optimization as a roots problem, a variety of direct numerical optimization methods are available. These methods are available for both one-dimensional and multidimensional problems. As the name implies, one-dimensional problems involve functions that depend on a single dependent variable. As in Fig. 7.3a,the search then consists of climbing or descending one-dimensional peaks and valleys. Multi- dimensional problems involve functions that depend on two or more dependent variables. FIGURE 7.3 (a) One-dimensional optimization. This figure also illustrates how minimization of f (x) is equivalent to the maximization of − f (x). (b) Two-dimensional optimization. Note that this figure can be taken to represent either a maximization (contours increase in elevation up to the maximum like a mountain) or a minimization (contours decrease in elevation down to the minimum like a valley). f (x) y Optimum f (x , y ) f (x) f (x, y) x Minimum f (x) y x Maximum – f (x) – f (x) x x (a) (b)cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 186 186 OPTIMIZATION In the same spirit, a two-dimensional optimization can again be visualized as searching out peaks and valleys (Fig. 7.3b). However, just as in real hiking, we are not constrained to walk a single direction; instead the topography is examined to efficiently reach the goal. Finally, the process of finding a maximum versus finding a minimum is essentially ∗ identical because the same value x both minimizes f (x) and maximizes − f (x). This equivalence is illustrated graphically for a one-dimensional function in Fig. 7.3a. In the next section, we will describe some of the more common approaches for one- dimensional optimization. Then we will provide a brief description of how MATLAB can be employed to determine optima for multidimensional functions. 7.2 ONE-DIMENSIONAL OPTIMIZATION This section will describe techniques to find the minimum or maximum of a function of a single variable f (x). A useful image in this regard is the one-dimensional “roller coaster”–like function depicted in Fig. 7.4. Recall from Chaps. 5 and 6 that root location was complicated by the fact that several roots can occur for a single function. Similarly, both local and global optima can occur in optimization. A global optimum represents the very best solution. A local optimum, though not the very best, is better than its immediate neighbors. Cases that include local optima are called multimodal. In such cases, we will almost always be interested in finding the global optimum. In addition, we must be concerned about mistaking a local result for the global optimum. Just as in root location, optimization in one dimension can be divided into bracketing and open methods. As described in the next section, the golden-section search is an example of a bracketing method that is very similar in spirit to the bisection method for root location. This is followed by a somewhat more sophisticated bracketing approach—parabolic inter- polation. We will then show how these two methods are combined and implemented with MATLAB’s fminbnd function. FIGURE 7.4 A function that asymptotically approaches zero at plus and minus ∞ and has two maximum and two minimum points in the vicinity of the origin. The two points to the right are local optima, whereas the two to the left are global. f(x) Local maximum Global maximum x Local Global minimum minimumcha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 187 7.2 ONE-DIMENSIONAL OPTIMIZATION 187 7.2.1 Golden-Section Search In many cultures, certain numbers are ascribed magical qualities. For example, we in the West are all familiar with “lucky 7” and “Friday the 13th.” Beyond such superstitious quantities, there are several well-known numbers that have such interesting and powerful mathematical properties that they could truly be called “magical.” The most common of these are the ratio of a circle’s circumference to its diameter π and the base of the natural logarithm e. Although not as widely known, the golden ratio should surely be included in the pan- theon of remarkable numbers. This quantity, which is typically represented by the Greek letter φ (pronounced: fee), was originally defined by Euclid (ca. 300 BCE) because of its role in the construction of the pentagram or five-pointed star. As depicted in Fig. 7.5, Euclid’s definition reads: “A straight line is said to have been cut in extreme and mean ratio when, as the whole line is to the greater segment, so is the greater to the lesser.” The actual value of the golden ratio can be derived by expressing Euclid’s definition as  +   1 2 1 = (7.2)   1 2 Multiplying by  / and collecting terms yields 1 2 2 φ − φ − 1 = 0 (7.3) where φ =  / . The positive root of this equation is the golden ratio: 1 2 √ 1 + 5 φ = = 1.61803398874989... (7.4) 2 The golden ratio has long been considered aesthetically pleasing in Western cultures. In addition, it arises in a variety of other contexts including biology. For our purposes, it provides the basis for the golden-section search, a simple, general-purpose method for de- termining the optimum of a single-variable function. The golden-section search is similar in spirit to the bisection approach for locating roots in Chap. 5. Recall that bisection hinged on defining an interval, specified by a lower guess (x ) and an upper guess (x ) that bracketed a single root. The presence of a root be- l u tween these bounds was verified by determining that f (x ) and f (x ) had different signs. l u The root was then estimated as the midpoint of this interval: x + x l u x = (7.5) r 2 FIGURE 7.5 Euclid’s definition of the golden ratio is based on dividing a line into two segments so that the ratio of the whole line to the larger segment is equal to the ratio of the larger segment to the smaller segment. This ratio is called the golden ratio.   1 2   1 2cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 188 188 OPTIMIZATION The final step in a bisection iteration involved determining a new smaller bracket. This was done by replacing whichever of the bounds x or x had a function value with the same sign l u as f (x ). A key advantage of this approach was that the new value x replaced one of the old r r bounds. Now suppose that instead of a root, we were interested in determining the minimum of a one-dimensional function. As with bisection, we can start by defining an interval that contains a single answer. That is, the interval should contain a single minimum, and hence is called unimodal. We can adopt the same nomenclature as for bisection, where x and x l u defined the lower and upper bounds, respectively, of such an interval. However, in contrast to bisection, we need a new strategy for finding a minimum within the interval. Rather than using a single intermediate value (which is sufficient to detect a sign change, and hence a zero), we would need two intermediate function values to detect whether a minimum occurred. The key to making this approach efficient is the wise choice of the intermediate points. As in bisection, the goal is to minimize function evaluations by replacing old values with new values. For bisection, this was accomplished by choosing the midpoint. For the golden-section search, the two intermediate points are chosen according to the golden ratio: x = x + d (7.6) 1 l x = x − d (7.7) 2 u where d = (φ − 1)(x − x ) (7.8) u l The function is evaluated at these two interior points. Two results can occur: 1. If, as in Fig. 7.6a, f (x ) f (x ), then f (x ) is the minimum, and the domain of x to the 1 2 1 left of x , from x to x , can be eliminated because it does not contain the minimum. For 2 l 2 this case, x becomes the new x for the next round. 2 l 2. If f (x ) f (x ), then f (x ) is the minimum and the domain of x to the right of x , from 2 1 2 1 x to x would be eliminated. For this case, x becomes the new x for the next round. 1 u 1 u Now, here is the real benefit from the use of the golden ratio. Because the original x 1 and x were chosen using the golden ratio, we do not have to recalculate all the function 2 values for the next iteration. For example, for the case illustrated in Fig. 7.6, the old x be- 1 comes the new x . This means that we already have the value for the new f (x ), since it is 2 2 the same as the function value at the old x . 1 To complete the algorithm, we need only determine the new x . This is done with 1 Eq. (7.6) with d computed with Eq. (7.8) based on the new values of x and x . A similar l u approach would be used for the alternate case where the optimum fell in the left subinterval. For this case, the new x would be computed with Eq. (7.7). 2 As the iterations are repeated, the interval containing the extremum is reduced rapidly. In fact, each round the interval is reduced by a factor of φ − 1 (about 61.8%). That means 10 that after 10 rounds, the interval is shrunk to about 0.618 or 0.008 or 0.8% of its initial length. After 20 rounds, it is about 0.0066%. This is not quite as good as the reduction achieved with bisection (50%), but this is a harder problem.cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 189 7.2 ONE-DIMENSIONAL OPTIMIZATION 189 f(x) Eliminate Minimum x x x l d 1 x x d 2 u (a) f(x) Old Old x x 2 1 x x x x x l 2 1 u New (b) FIGURE 7.6 (a) The initial step of the golden-section search algorithm involves choosing two interior points according to the golden ratio. (b) The second step involves defining a new interval that encompasses the optimum. EXAMPLE 7.2 Golden-Section Search Problem Statement. Use the golden-section search to find the minimum of 2 x f (x) = − 2 sin x 10 within the interval from x = 0 to x = 4. l u Solution. First, the golden ratio is used to create the two interior points: d = 0.61803(4 − 0) = 2.4721 x = 0 + 2.4721 = 2.4721 1 x = 4 − 2.4721 = 1.5279 2cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 190 190 OPTIMIZATION The function can be evaluated at the interior points: 2 1.5279 f (x ) = − 2 sin(1.5279)=−1.7647 2 10 2 2.4721 f (x ) = − 2 sin(2.4721)=−0.6300 1 10 Because f (x ) f (x ), our best estimate of the minimum at this point is that it is 2 1 located at x = 1.5279 with a value of f (x) = –1.7647. In addition, we also know that the minimum is in the interval defined by x , x , and x . Thus, for the next iteration, the lower l 2 1 bound remains x = 0, and x becomes the upper bound, that is, x = 2.4721. In addition, l 1 u the former x value becomes the new x , that is, x = 1.5279. In addition, we do not have to 2 1 1 recalculate f (x ), it was determined on the previous iteration as f (1.5279) = –1.7647. 1 All that remains is to use Eqs. (7.8) and (7.7) to compute the new value of d and x : 2 d = 0.61803(2.4721 − 0) = 1.5279 x = 2.4721 − 1.5279 = 0.9443 2 The function evaluation at x is f (0.9943)=−1.5310. Since this value is less than the 2 function value at x , the minimum is f (1.5279)=−1.7647, and it is in the interval pre- 1 scribed by x , x , and x . The process can be repeated, with the results tabulated here: 2 1 u ix f (x ) x f (x ) x f (x ) x f (x ) d l l 2 2 1 1 u u 1 0 0 1.5279 −1.7647 2.4721 −0.6300 4.0000 3.1136 2.4721 2 0 0 0.9443 −1.5310 1.5279 −1.7647 2.4721 −0.6300 1.5279 3 0.9443 −1.5310 1.5279 −1.7647 1.8885 −1.5432 2.4721 −0.6300 0.9443 4 0.9443 −1.5310 1.3050 −1.7595 1.5279 −1.7647 1.8885 −1.5432 0.5836 5 1.3050 −1.7595 1.5279 −1.7647 1.6656 −1.7136 1.8885 −1.5432 0.3607 6 1.3050 −1.7595 1.4427 −1.7755 1.5279 −1.7647 1.6656 −1.7136 0.2229 7 1.3050 −1.7595 1.3901 −1.7742 1.4427 −1.7755 1.5279 −1.7647 0.1378 8 1.3901 −1.7742 1.4427 −1.7755 1.4752 −1.7732 1.5279 −1.7647 0.0851 Note that the current minimum is highlighted for every iteration. After the eighth iteration, the minimum occurs at x = 1.4427 with a function value of −1.7755. Thus, the result is converging on the true value of −1.7757 at x = 1.4276. Recall that for bisection (Sec. 5.4), an exact upper bound for the error can be calcu- lated at each iteration. Using similar reasoning, an upper bound for golden-section search can be derived as follows: Once an iteration is complete, the optimum will either fall in one of two intervals. If the optimum function value is at x , it will be in the lower interval (x , 2 l x , x ). If the optimum function value is at x , it will be in the upper interval (x , x , x ). 2 1 1 2 1 u Because the interior points are symmetrical, either case can be used to define the error.cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 191 7.2 ONE-DIMENSIONAL OPTIMIZATION 191 Looking at the upper interval (x , x , x ), if the true value were at the far left, the max- 2 1 u imum distance from the estimate would be x = x − x a 1 2 = x + (φ − 1)(x − x ) − x + (φ − 1)(x − x ) l u l u u l = (x − x ) + 2(φ − 1)(x − x ) l u u l = (2φ − 3)(x − x ) u l or 0.2361 (x − x ). If the true value were at the far right, the maximum distance from the u l estimate would be x = x − x b u 1 = x − x − (φ − 1)(x − x ) u l u l = (x − x ) − (φ − 1)(x − x ) u l u l = (2 − φ)(x − x ) u l or 0.3820 (x − x ). Therefore, this case would represent the maximum error. This result can u l then be normalized to the optimal value for that iteration x to yield opt     x − x u l   ε = (2 − φ) × 100% (7.9) a   x opt This estimate provides a basis for terminating the iterations. An M-file function for the golden-section search for minimization is presented in Fig. 7.7. The function returns the location of the minimum, the value of the function, the approximate error, and the number of iterations. The M-file can be used to solve the problem from Example 7.1. g=9.81;v0=55;m=80;c=15;z0=100; z=(t) -(z0+m/c(v0+mg/c)(1-exp(-c/mt))-mg/ct); xmin,fmin,ea,iter=goldmin(z,0,8) xmin = 3.8317 fmin = -192.8609 ea = 6.9356e-005 Notice how because this is a maximization, we have entered the negative of Eq. (7.1). Consequently, fmin corresponds to a maximum height of 192.8609. You may be wondering why we have stressed the reduced function evaluations of the golden-section search. Of course, for solving a single optimization, the speed savings would be negligible. However, there are two important contexts where minimizing the number of function evaluations can be important. These are 1. Many evaluations. There are cases where the golden-section search algorithm may be a part of a much larger calculation. In such cases, it may be called many times. Therefore, keeping function evaluations to a minimum could pay great dividends for such cases.cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 192 192 OPTIMIZATION function x,fx,ea,iter=goldmin(f,xl,xu,es,maxit,varargin) % goldmin: minimization golden section search % x,fx,ea,iter=goldmin(f,xl,xu,es,maxit,p1,p2,...): % uses golden section search to find the minimum of f % input: % f = name of function % xl, xu = lower and upper guesses % es = desired relative error (default = 0.0001%) % maxit = maximum allowable iterations (default = 50) % p1,p2,... = additional parameters used by f % output: % x = location of minimum % fx = minimum function value % ea = approximate relative error (%) % iter = number of iterations if nargin3,error('at least 3 input arguments required'),end if nargin4isempty(es), es=0.0001;end if nargin5isempty(maxit), maxit=50;end phi=(1+sqrt(5))/2; iter=0; while(1) d = (phi-1)(xu - xl); x1 = xl + d; x2 = xu - d; if f(x1,varargin:) f(x2,varargin:) xopt = x1; xl = x2; else xopt = x2; xu = x1; end iter=iter+1; if xopt=0, ea = (2 - phi) abs((xu - xl) / xopt) 100;end if ea = es iter = maxit,break,end end x=xopt;fx=f(xopt,varargin:); FIGURE 7.7 An M-file to determine the minimum of a function with the golden-section search. 2. Time-consuming evaluation. For pedagogical reasons, we use simple functions in most of our examples. You should understand that a function can be very complex and time-consuming to evaluate. For example, optimization can be used to estimate the parameters of a model consisting of a system of differential equations. For such cases, the “function” involves time-consuming model integration. Any method that minimizes such evaluations would be advantageous.cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 193 7.2 ONE-DIMENSIONAL OPTIMIZATION 193 Parabolic True maximum approximation f(x) of maximum True function Parabolic function x x x x x 1 2 4 3 FIGURE 7.8 Graphical depiction of parabolic interpolation. 7.2.2 Parabolic Interpolation Parabolic interpolation takes advantage of the fact that a second-order polynomial often provides a good approximation to the shape of f(x) near an optimum (Fig. 7.8). Just as there is only one straight line connecting two points, there is only one parabola connecting three points. Thus, if we have three points that jointly bracket an optimum, we can fit a parabola to the points. Then we can differentiate it, set the result equal to zero, and solve for an estimate of the optimal x. It can be shown through some algebraic manipula- tions that the result is 2 2 1 (x − x ) f (x ) − f (x ) − (x − x ) f (x ) − f (x ) 2 1 2 3 2 3 2 1 x = x − (7.10) 4 2 2 (x − x ) f (x ) − f (x ) − (x − x ) f (x ) − f (x ) 2 1 2 3 2 3 2 1 where x , x , and x are the initial guesses, and x is the value of x that corresponds to the 1 2 3 4 optimum value of the parabolic fit to the guesses. EXAMPLE 7.3 Parabolic Interpolation Problem Statement. Use parabolic interpolation to approximate the minimum of 2 x f (x) = − 2 sin x 10 with initial guesses of x = 0, x = 1, and x = 4. 1 2 3 Solution. The function values at the three guesses can be evaluated: x = 0 f(x ) = 0 1 1 x = 1 f(x ) =−1.5829 2 2 x = 4 f(x ) = 3.1136 3 3cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 194 194 OPTIMIZATION and substituted into Eq. (7.10) to give 2 2 1 (1 − 0) −1.5829 − 3.1136 − (1 − 4) −1.5829 − 0 x = 1 − = 1.5055 4 2 (1 − 0) −1.5829 − 3.1136 − (1 − 4) −1.5829 − 0 which has a function value of f(1.5055) =−1.7691. Next, a strategy similar to the golden-section search can be employed to determine which point should be discarded. Because the function value for the new point is lower than for the intermediate point (x ) and the new x value is to the right of the intermediate 2 point, the lower guess (x ) is discarded. Therefore, for the next iteration: 1 x = 1 f(x ) =−1.5829 1 1 x = 1.5055 f(x ) =−1.7691 2 2 x = 4 f(x ) = 3.1136 3 3 which can be substituted into Eq. (7.10) to give 2 2 1 (1.5055 − 1) −1.7691 − 3.1136 − (1.5055 − 4) −1.7691 − (−1.5829) x = 1.5055 − 4 2 (1.5055 − 1) −1.7691 − 3.1136 − (1.5055 − 4) −1.7691 − (−1.5829) = 1.4903 which has a function value of f(1.4903) =−1.7714. The process can be repeated, with the results tabulated here: ix f (x ) x f (x ) x f (x ) x f (x ) 1 1 2 2 3 3 4 4 1 0.0000 0.0000 1.0000 −1.5829 4.0000 3.1136 1.5055 −1.7691 2 1.0000 −1.5829 1.5055 −1.7691 4.0000 3.1136 1.4903 −1.7714 3 1.0000 −1.5829 1.4903 −1.7714 1.5055 −1.7691 1.4256 −1.7757 4 1.0000 −1.5829 1.4256 −1.7757 1.4903 −1.7714 1.4266 −1.7757 5 1.4256 −1.7757 1.4266 −1.7757 1.4903 −1.7714 1.4275 −1.7757 Thus, within five iterations, the result is converging rapidly on the true value of −1.7757 at x = 1.4276. 7.2.3 MATLAB Function: fminbnd Recall that in Sec. 6.4 we described Brent’s method for root location, which combined sev- eral root-finding methods into a single algorithm that balanced reliability with efficiency. Because of these qualities, it forms the basis for the built-in MATLAB function fzero. Brent also developed a similar approach for one-dimensional minimization which forms the basis for the MATLAB fminbnd function. It combines the slow, dependable golden-section search with the faster, but possibly unreliable, parabolic interpolation. It first attempts parabolic interpolation and keeps applying it as long as acceptable results are obtained. If not, it uses the golden-section search to get matters in hand.cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 195 7.3 MULTIDIMENSIONAL OPTIMIZATION 195 A simple expression of its syntax is xmin, fval = fminbnd(function,x1,x2) where x and fval are the location and value of the minimum, function is the name of the function being evaluated, and x1 and x2 are the bounds of the interval being searched. Here is a simple MATLAB session that uses fminbnd to solve the problem from Example 7.1. g=9.81;v0=55;m=80;c=15;z0=100; z=(t) -(z0+m/c(v0+mg/c)(1-exp(-c/mt))-mg/ct); x,f=fminbnd(z,0,8) x = 3.8317 f = -192.8609 As with fzero, optional parameters can be specified using optimset. For example, we can display calculation details: options = optimset('display','iter'); fminbnd(z,0,8,options) Func-count x f(x) Procedure 1 3.05573 -189.759 initial 2 4.94427 -187.19 golden 3 1.88854 -171.871 golden 4 3.87544 -192.851 parabolic 5 3.85836 -192.857 parabolic 6 3.83332 -192.861 parabolic 7 3.83162 -192.861 parabolic 8 3.83166 -192.861 parabolic 9 3.83169 -192.861 parabolic Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 ans = 3.8317 Thus, after three iterations, the method switches from golden to parabolic, and after eight iterations, the minimum is determined to a tolerance of 0.0001. 7.3 MULTIDIMENSIONAL OPTIMIZATION Aside from one-dimensional functions, optimization also deals with multidimensional functions. Recall from Fig. 7.3a that our visual image of a one-dimensional search was like a roller coaster. For two-dimensional cases, the image becomes that of mountains and valleys (Fig. 7.3b). As in the following example, MATLAB’s graphic capabilities provide a handy means to visualize such functions.cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 196 196 OPTIMIZATION EXAMPLE 7.4 Visualizing a Two-Dimensional Function Problem Statement. Use MATLAB’s graphical capabilities to display the following function and visually estimate its minimum in the range –2 ≤ x ≤ 0 and 0 ≤ x ≤ 3: 1 2 2 2 f (x , x ) = 2 + x − x + 2x + 2x x + x 1 2 1 2 1 2 1 2 Solution. The following script generates contour and mesh plots of the function: x=linspace(-2,0,40);y=linspace(0,3,40); X,Y = meshgrid(x,y); Z=2+X-Y+2X.2+2X.Y+Y.2; subplot(1,2,1); cs=contour(X,Y,Z);clabel(cs); xlabel('x_1');ylabel('x_2'); title('(a) Contour plot');grid; subplot(1,2,2); cs=surfc(X,Y,Z); zmin=floor(min(Z)); zmax=ceil(max(Z)); xlabel('x_1');ylabel('x_2');zlabel('f(x_1,x_2)'); title('(b) Mesh plot'); As displayed in Fig. 7.9, both plots indicate that function has a minimum value of about f(x , x ) = 0 to 1 located at about x =−1 and x = 1.5. 1 2 1 2 FIGURE 7.9 (a) Contour and (b) mesh plots of a two-dimensional function. (a) Contour plot (b) Mesh plot 3 7  3 6   2.5 8 4   1  6 2 2  4 1.5 2 3 1  0 3 0 0.5 2 2  4  –1 x 1 2 x 7 5 8 6 1    0 0 –2 –2 –1.5 –1 –0.5 0 x 1 x 2 f (x , x ) 1 2cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 197 7.4 CASE STUDY 197 Techniques for multidimensional unconstrained optimization can be classified in a number of ways. For purposes of the present discussion, we will divide them depending on whether they require derivative evaluation. Those that require derivatives are called gradient, or descent (or ascent), methods. The approaches that do not require derivative eval- uation are called nongradient, or direct, methods. As described next, the built-in MATLAB functionfminsearch is a direct method. 7.3.1 MATLAB Function: fminsearch Standard MATLAB has a function fminsearch that can be used to determine the mini- mum of a multidimensional function. It is based on the Nelder-Mead method, which is a direct-search method that uses only function values (does not require derivatives) and handles non-smooth objective functions. A simple expression of its syntax is xmin, fval = fminsearch(function,x0) wherexmin andfval are the location and value of the minimum,function is the name of the function being evaluated, andx0 is the initial guess. Note that x0 can be a scalar, vector, or a matrix. Here is a simple MATLAB session that uses fminsearch to determine minimum for the function we just graphed in Example 7.4: f=(x) 2+x(1)-x(2)+2x(1)2+2x(1)x(2)+x(2)2; x,fval=fminsearch(f,-0.5,0.5) x = -1.0000 1.5000 fval = 0.7500 7.4 CASE STUDY EQUILIBRIUM AND MINIMUM POTENTIAL ENERGY Background. As in Fig. 7.10a, an unloaded spring can be attached to a wall mount. When a horizontal force is applied, the spring stretches. The displacement is related to the force by Hookes law, F = kx. The potential energy of the deformed state consists of the dif- ference between the strain energy of the spring and the work done by the force: (7.11) 2 PE(x) = 0.5kx − Fx FIGURE 7.10 (a) An unloaded spring attached to a wall mount. (b) Application of a horizontal force stretches the spring where the relationship between force and displacement is described by Hooke’s law. k (a) x F (b)cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 198 198 OPTIMIZATION 7.4 CASE STUDY continued Equation (7.11) defines a parabola. Since the potential energy will be at a minimum at equilibrium, the solution for displacement can be viewed as a one-dimensional optimiza- tion problem. Because this equation is so easy to differentiate, we can solve for the dis- placement as x = F/k. For example, if k = 2 N/cm and F = 5 N, x = 5N/(2 N/cm) = 2.5 cm. A more interesting two-dimensional case is shown in Fig. 7.11. In this system, there are two degrees of freedom in that the system can move both horizontally and vertically. In the same way that we approached the one-dimensional system, the equilibrium defor- mations are the values of x and x that minimize the potential energy: 1 2   2 2 2 PE(x , x ) = 0.5k x + (L − x ) − L 1 2 a a 2 a 1   2 2 2 + 0.5k x + (L + x ) − L − F x − F x (7.12) b b 2 b 1 1 2 2 1 If the parameters are k = 9 N/cm, k = 2 N/cm, L = 10 cm, L = 10 cm, F = 2 N, and a b a b 1 F = 4 N, use MATLAB to solve for the displacements and the potential energy. 2 FIGURE 7.11 A two-spring system: (a) unloaded and (b) loaded. F 2 L k a a F 1 x 2 x 1 L k b b (a) (b)cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 199 PROBLEMS 199 7.4 CASE STUDY continued Solution. An M-file can be developed to hold the potential energy function: function p=PE(x,ka,kb,La,Lb,F1,F2) PEa=0.5ka(sqrt(x(1)2+(La-x(2))2)-La)2; PEb=0.5kb(sqrt(x(1)2+(Lb+x(2))2)-Lb)2; W=F1x(1)+F2x(2); p=PEa+PEb-W; The solution can be obtained with the fminsearch function: ka=9;kb=2;La=10;Lb=10;F1=2;F2=4; x,f=fminsearch(PE,-0.5,0.5,,ka,kb,La,Lb,F1,F2) x = 4.9523 1.2769 f = -9.6422 Thus, at equilibrium, the potential energy is −9.6422 N·cm. The connecting point is located 4.9523 cm to the right and 1.2759 cm above its original position. PROBLEMS 7.1 Perform three iterations of the Newton-Raphson method (c) Differentiate the function and then use a root-location to determine the root of Eq. (E7.1.1). Use the parameter val- method to solve for the maximum f(x) and the corre- ues from Example 7.1 along with an initial guess of t = 3 s. sponding value of x. 7.2 Given the formula 7.5 Solve for the value of x that maximizes f(x) in Prob. 7.4 using the golden-section search. Employ initial guesses of 2 f (x)=−x + 8x − 12 x = 0 and x = 2, and perform three iterations. l u (a) Determine the maximum and the corresponding value of 7.6 Repeat Prob. 7.5, except use parabolic interpolation. x for this function analytically (i.e., using differentiation). Employ initial guesses of x = 0, x = 1, and x = 2, and per- 1 2 3 (b) Verify that Eq. (7.10) yields the same results based on form three iterations. initial guesses of x = 0, x = 2, and x = 6. 7.7 Employ the following methods to find the maximum of 1 2 3 7.3 Consider the following function: 2 3 4 f (x) = 4x − 1.8x + 1.2x − 0.3x 2 3 4 f (x) = 3 + 6x + 5x + 3x + 4x (a) Golden-section search (x = – 2, x = 4, ε = 1%). l u s Locate the minimum by finding the root of the derivative of (b) Parabolic interpolation (x = 1.75, x = 2, x = 2.5, 1 2 3 this function. Use bisection with initial guesses of x =−2 l iterations = 5). and x = 1. u 7.8 Consider the following function: 7.4 Given 4 3 2 6 4 f (x) = x + 2x + 8x + 5x f (x)=−1.5x − 2x + 12x (a) Plot the function. Use analytical and graphical methods to show the function (b) Use analytical methods to prove that the function is con- has a minimum for some value of x in the range cave for all values of x. −2 ≤ x ≤ 1.cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 200 200 OPTIMIZATION 7.9 Employ the following methods to find the minimum of T air the function from Prob. 7.8: r i (a) Golden-section search (x =−2, x = 1, ε = 1%). l u s (b) Parabolic interpolation (x =−2, x =−1, x = 1, 1 2 3 r w iterations = 5). 7.10 Consider the following function: T w 3 f (x) = 2x + x Perform 10 iterations of parabolic interpolation to locate the minimum. Comment on the convergence of your results (x = 0.1, x = 0.5, x = 5) 1 2 3 7.11 Develop a single script to (a) generate contour and FIGURE P7.14 mesh subplots of the following temperature field in a similar Cross-section of an insulated wire. fashion to Example 7.4: 2 2 T(x, y) = 2x + 3y − 4xy − y − 3x and (b) determine the minimum with fminsearch. in testing the response of the truss to a force exerted in any 7.12 The head of a groundwater aquifer is described in number of directions designated by the angle θ. The Cartesian coordinates by parameters for the problem are E = Young’s modulus 1 2 11 = 2 × 10 Pa, A = cross-sectional area = 0.0001 m , w = h(x, y) = 2 2 1 + x + y + x + xy width = 0.44 m,  = length = 0.56 m, and h = height = 0.5 m. The displacements x and y can be solved by deter- Develop a single script to (a) generate contour and mesh mining the values that yield a minimum potential energy. subplots of the function in a similar fashion to Example 7.4, Determine the displacements for a force of 10,000 N and a and (b) determine the maximum with fminsearch. range of θ’s from 0° (horizontal) to 90° (vertical). 7.13 Recent interest in competitive and recreational cycling 7.14 As electric current moves through a wire (Fig. P7.14), has meant that engineers have directed their skills toward the heat generated by resistance is conducted through a layer of design and testing of mountain bikes (Fig. P7.13a). Suppose insulation and then convected to the surrounding air. The that you are given the task of predicting the horizontal steady-state temperature of the wire can be computed as and vertical displacement of a bike bracketing system in   response to a force. Assume the forces you must analyze can q 1 r + r 1 1 w i be simplified as depicted in Fig. P7.13b. You are interested T = T + ln + air 2π k r h r + r w w i Determine the thickness of insulation r (m) that minimizes i the wire’s temperature given the following parameters: q = heat generation rate = 75 W/m, r = wire radius = 6 mm, w w k = thermal conductivity of insulation = 0.17 W/(m K), 2 h = convective heat transfer coefficient = 12 W/(m K), and T = air temperature = 293 K. air 7.15 Develop an M-file that is expressly designed to locate h  a maximum with the golden-section search. In other words, x set if up so that it directly finds the maximum rather than  finding the minimum of −f(x). The function should have the F y following features: (a) (b) • Iterate until the relative error falls below a stopping cri- terion or exceeds a maximum number of iterations. FIGURE P7.13 • Return both the optimal x and f(x). (a) A mountain bike along with (b) a free-body diagram for a part of the frame. Test your program with the same problem as Example 7.1.cha01102_ch07_182-204.qxd 12/17/10 8:06 AM Page 201 PROBLEMS 201 7.22 The normal distribution is a bell-shaped curve defined by 7.16 Develop an M-file to locate a minimum with the 2 golden-section search. Rather than using the maximum itera- −x y = e tions and Eq. (7.9) as the stopping criteria, determine the Use the golden-section search to determine the location of number of iterations needed to attain a desired tolerance. Test the inflection point of this curve for positive x. your function by solving Example 7.2 using E = 0.0001. a,d 7.23 Use the fminsearch function to determine the 7.17 Develop an M-file to implement parabolic interpola- minimum of tion to locate a minimum. The function should have the following features: 2 2 f (x, y) = 2y − 2.25xy − 1.75y + 1.5x • Base it on two initial guesses, and have the program 7.24 Use the fminsearch function to determine the generate the third initial value at the midpoint of the maximum of interval. 2 4 2 • Check whether the guesses bracket a maximum. If not, f (x, y) = 4x + 2y + x − 2x + 2xy − 3y the function should not implement the algorithm, but 7.25 Given the following function: should return an error message. • Iterate until the relative error falls below a stopping cri- 2 2 f (x, y)=−8x + x + 12y + 4y − 2xy terion or exceeds a maximum number of iterations. • Return both the optimal x and f(x). Determine the minimum (a) graphically, (b) numerically with the fminsearch function, and (c) substitute the result Test your program with the same problem as Example 7.3. of (b) back into the function to determine the minimum 7.18 Pressure measurements are taken at certain points f (x, y). behind an airfoil over time. These data best fit the curve 7.26 The specific growth rate of a yeast that produces an y = 6 cos x − 1.5 sin x from x = 0 to 6 s. Use four iterations antibiotic is a function of the food concentration c: of the golden-search method to find the minimum pressure. Set x = 2 and x = 4. 2c l u g = 2 3 7.19 The trajectory of a ball can be computed with 4 + 0.8c + c + 0.2c g 2 y = (tanθ )x − x + y As depicted in Fig. P7.26, growth goes to zero at very low 0 0 2 2 2v cos θ 0 0 concentrations due to food limitation. It also goes to zero at high concentrations due to toxicity effects. Find the value of where y = the height (m), θ = the initial angle (radians), 0 c at which growth is a maximum. v = the initial velocity (m/s), g = the gravitational 0 2 7.27 A compound A will be converted into B in a stirred constant = 9.81 m/s , and y = the initial height (m). Use the 0 tank reactor. The product B and unreacted A are purified in a golden-section search to determine the maximum height ◦ separation unit. Unreacted A is recycled to the reactor. A given y = 1 m, v = 25 m/s, and θ = 50 . Iterate until the 0 0 0 process engineer has found that the initial cost of the system approximate error falls below ε = 1% using initial guesses s of x = 0 and x = 60 m. l u 7.20 The deflection of a uniform beam subject to a linearly increasing distributed load can be computed as 0.4 w 0 5 2 3 4 y = (−x + 2L x − L x) 120EIL g 1 0.2 (d ) 2 4 Given that L = 600 cm, E = 50,000 kN/cm , I = 30,000 cm , and w = 2.5 kN/cm, determine the point of maximum de- 0 flection (a) graphically, (b) using the golden-section search 0 until the approximate error falls below ε = 1% with initial 0 10 5 s guesses of x = 0 and x = L. c (mg/L) l u 7.21 A object with a mass of 90 kg is projected upward from the surface of the earth at a velocity of 60 m/s. If the FIGURE P7.26 object is subject to linear drag (c = 15 kg/s), use the golden- The specific growth rate of a yeast that produces an section search to determine the maximum height the object antibiotic versus the food concentration. attains.

Advise: Why You Wasting Money in Costly SEO Tools, Use World's Best Free SEO Tool Ubersuggest.