m

All materials on our website are shared by users. If you have any questions about copyright issues, please report us to resolve them. We are always happy to assist you.

Description

CEE 201L. Uncertainty, Design, and Optimization Department of Civil and Environmental Engineering Duke University Homework 3, due: Friday January 31, 2020 In this course, we will apply numerical methods

Transcript

CEE 201L. Uncertainty, Design, and Optimization Department of Civil and Environmental Engineering Duke University Homework 3, due: Friday January 31, 2020 In this course, we will apply numerical methods to solve constrained optimization problems. Mathematically speaking, in a constrained optimization problem we are asked to find values for n [ ] T, parameters or design variables, x = x 1 x 2 x 3 x n that minimizes a specified objective function f(x) such that a set of inequalities are satisfied. In many design problems, the objective is to simply minimize the total cost (or price) of the product you are designing. In other design problems, the objective could be to maximize an aspect of the performance, to maximize profits, etc. The suitability of the optimized design depends on more than the cost (or performance, profits, etc). Most systems must also meet a range of other criteria (e.g., safe enough, reliable enough, strong enough). These criteria are expressed as inequality constraints, and depend upon the set of design variables. As long as the inequalities are met the criteria are satisfied. The inequalities constrain the best value of the objective. In a structural design, for example, some inequality constraints might be that: (i) the stress remains less than a yield stress value, (ii) the compression force in a column remains less than the buckling load of the column, and (iii) the deformation of the system remains within certain limits. Note that any inequality a b can be written a b 0 or a/b 1 0. In general, we will write a set of m inequality constraints as g(x) 0, or g 1 (x) 0, g 2 (x) 0,, g m (x) 0. With this notation, a constrained optimization problem may be written: minimize x 1, x 2,..., x n f(x 1, x 2,..., x n ), such that g 1 (x 1, x 2,..., x n ) 0 g 2 (x 1, x 2,..., x n ) 0. g m (x 1, x 2,..., x n ) 0 (1) In certain rare cases, we can write actual equations for f(x) and g(x) 0, and use calculus to find an equation for the constrained optimum. In the vast majority of practical problems, however, the equations are simply much too complicated for this approach. In such cases computer-aided analysis can automate the evaluation of the cost and feasibility of a particular trial design. Further, computer-aided optimization allows designers to rapidly determine feasible and possibly optimal designs. In this homework, you will examine this process by solving optimization problems using a number of methods and will investigate the mathematics of quadratic programming. In the first problem, you will see that optimization problems with discontinuous feasible regions are not hard to come up with. The global minimum of such problems can be found with a gridded search or by repeatedly using local minimization methods from different initial guesses. In the second problem, you will see that linearly-constrained quadratic objectives can be optimized by solving a simple linear matrix equation. For such problems, no initial guess is required, and if there is a feasible solution, there is only one feasible region. This problem will also provide some insight that can be useful for interpreting the values of the Lagrange multipliers. 2 Homework 3 A gridded search in 2D A continuous parameter space (x 1, x 2 ) can be sampled into a finite number of points (shown as blue circles in the figure below) on a grid from x 1,min to x 1,max and from x 2,min to x 2,max. The cost f(x 1, x 2 ) and constraints g(x 1, x 2 ) can be evaluated at every point on the grid. By keeping track of the best feasible point as each point is evaluated, the optimal solution from within the set of grid points can be identified (the green star). In most cases, the true optimal value (black star) is off of the grid. A finer grid would contain a point closer to the true optimal, but would also involve more evaluations of f and g. This idea raises the idea of a method to adaptively refine the grid around (one or more) approximate optimal solutions. The figure below is representative of problems that have discontinuous feasible regions. A local minimizing routine (such as SQP) starting within region A will never move into region B as this would require moving through the infeasible space (hashed). Search methods in which constraints are incorporated through a penalty method could potentially move from one feasible region into another, if the penalty factor is small enough. In 3D, the gridded space would be a cube divided into little blocks. The feasible regions would be 3D objects within the 3D grid. The grid shown below is a regular grid in which all values of x i are uniformly spaced with the same increment. x x 2,max 2,min A B x 1,min increment x 1,max CEE 201L. Uncertainty, Design, and Optimization Duke University 2020 H.P. Gavin 3 1. Constrained Optimization using a Gridded Parameter Search (25 points): This example optimization problem is to find the minimum of Fletcher and Powell s helical valley 1, which is a function of three parameters x 1, x 2, and x 3 { f(x 1, x 2, x 3 ) = 100 [x 3 10q(x 1, x 2 )] 2 + [r(x 1, x 2 ) 1] 2} + x 2 3 (2) where and and such that (2π)q(x 1, x 2 ) = arctan(x 2, x 1 ) (use atan2(x(2), x(1)); ) (3) r 2 (x 1, x 2 ) = x x 2 2 g 1 (x 1, x 2, x 3 ) = 3x 2 5x and (4) g 2 (x 1, x 2, x 3 ) = 10x 1 20 cos(2x 3 ) 0. (5) To further constrain the problem consider parameter values only within the bounds: 10 x 1 10, 10 x 2 10, and 2.5 x Write a Matlab function and a Matlab script. The Matlab function uses a trial set of values for x 1, x 2, and x 3 to calculate f(x 1, x 2, x 3 ), g 1 (x 1, x 2, x 3 ), and g 2 (x 1, x 2, x 3 ). This file name must be HW3_2020_analysis.m, and the first line of this function file must be: function [f,g] = HW3_2020_analysis(x,c); This function uses values in a three-element parameter vector x to compute a value of the cost function (the scalar variable f), and values of the constraints (the column-vector g). In this problem, the Matlab (column) vector g has two elements: [g(1) ; g(2)] This matlab function has just six lines of code. The Matlab script is a file containing Matlab commands that use the analysis function to find the optimal value of the parameters. This Matlab script, called HW3_2020_opt.m, is shown on the next page. It implements a method that is terribly inefficient, but is easy-to-understand, easy-to-program, and is guaranteed to find a value that is close to the optimal value. The script starts by defining a set of values for x 1, x 2, and x 3 that are within their limiting values (lines 1-8). It then determines the number of values in each set with the length command (lines 8-12). The total number of analyses is the product of the length of each parameter-set vector (line 13). Given sets of values for x1, x2, and x3, the script uses a set of nested for-loops to evaluate your function [f,g] = HW3_2020_analysis(x,c); for every combination of parameters (lines 22-49). For each evaluation, the script checks to see if all elements of g are less than zero and if f is less than the smallest value of f found so far (line 30). If all elements of g are less than zero (feasible) and f is less than the best value computed so far, the script sets f_opt = f; and x_opt = x; (lines 31-35). If not, the script simply goes on to the next combination of parameters. After all combinations have been tried, the optimal values for the parameters will be x_opt and the associated value of the cost function will be f_opt (lines 51-54). 1 Fletcher, R., and Powell, M.J.D. (1963). A rapidly convergent descent method for minimization, The Computer Journal (6):163. 4 Homework 3 Since this method evaluates lots and lots of alternative designs, you may expect the script to take some time to run. The script keeps track of how much time it takes using the commands tic and toc (lines 21 and 50). When you run the script Matlab will show how many seconds elapsed between tic and toc 1 x_min = [ ]; % minimum v a l u e s f o r each parameter 2 x_max = [ ]; % maximum v a l u e s f o r each parameter 3 4 increment = 1.0; % s m a l l e r increment, more v a l u e s to t r y 5 x1_set = [ x_min (1) : increment : x_max (1) ]; % t h e s e t o f v a l u e s f o r x (1) 6 x2_set = [ x_min (2) : increment : x_max (2) ]; % t h e s e t o f v a l u e s f o r x (2) 7 x3_set = [ x_min (3) : increment : x_max (3) ]; % t h e s e t o f v a l u e s f o r x (3) 8 9 N_x1 = length( x1_set ); 10 N_x2 = length( x2_set ); 11 N_x3 = length( x3_set ); NumberOfAnalyses = N_x1 * N_x2 * N_x f_opt = 1 e16 ; % i n i t i a l v a l u e f o r o b j e c t i v e f... some b i g number 16 g_opt = [1;1]; % i n i t i a l v a l u e s o f c o n s t r a i n t s more off ; % show r e s u l t s as t h e y are computed count = 0; 21 t i c 22 for i =1: N_x1 % l o o p over v a l u e s o f x (1) 23 for j =1: N_x2 % l o o p over v a l u e s o f x (2) 24 for k =1: N_x3 % l o o p over v a l u e s o f x (3) x = [ x1_set (i) ; x2_set (j) ; x3_set (k) ]; % t r i a l v a l u e s [f,g] = HW3_2020_analysis (x,1); % run t h e a n a l y s i s i f a l l (g 0) & f f_opt % b e s t s o l u t i o n, so f a r f_opt = f; % update o b j e c t i v e 33 x_opt = x; % update v a r i a b l e s 34 g_opt = g; % update c o n s t r a i n t s end count = count + 1; 39 i f rem( count,100) % when w i l l t h i s f i n i s h? 40 secs = toc; 41 secs_left = round(( NumberOfAnalyses - count )* secs / count ); 42 eta = datestr ( now + secs_left /3600/24,14); 43 f p r i n t f ( %4d (%4.1 f %%); %5.2 f secs / count ; eta : %s (%4 d s)\n, count,100* count / NumberOfAnalyses, secs / count, eta, secs_left ); 45 end end % x (3) l o o p 48 end % x (2) l o o p 49 end % x (1) l o o p 50 time_a = toc x_opt_a = x_opt % d i s p l a y t h e optimal parameters 53 f_opt_a = f_opt % d i s p l a y t h e optimal c o s t 54 g_opt_a = g_opt % d i s p l a y t h e c o n s t r a i n t s Consider everything. Keep the good. Avoid evil whenever you notice it. CEE 201L. Uncertainty, Design, and Optimization Duke University 2020 H.P. Gavin 5 To hand in: A print out of your well-commented Matlab function HW analysis.m Every Matlab file that you write should be well-commented. Immediately beneath the function [results] = my_func(variables) line write comments that tell a user about your m-file... what it does, how to use it, your name, your , and the date. The tables below, filled in. Do two analyses, first with increment = 1.0 and again with increment = 1.0. In each row, circle the value of the constraint(s) (g 1, g 2 ) that bind(s) the optimum point. Gridded Search Increment Value = 1.0 optimal values tic-toc f g 1 g 2 x 1 x 2 x 3 time Gridded Search Increment Value = 0.1 optimal values tic-toc f g 1 g 2 x 1 x 2 x 3 time 6 Homework 3 2. Mathematics of Quadratic Programming (30 points): This question involves matrix mathematics. The document linked here... hpgavin/cee201/matrixmath.pdf... provides some background on this topic. Consider a quadratic cost function in terms of three variable parameters x 1, x 2, x 3, f(x 1, x 2, x 3 ) = 52x x x x 1 x x 1 x 3 + 6x 2 x x x x (a) Derive three equations for f/ x 1, f/ x 2, and f/ x 3. (b) This cost function may be written in a form f(x) = 1 2 xt Hx + c T x + d (6) where H is a 3 3 Hessian matrix and c is a 3 1 vector and d is a constant. The Hessian matrix represents the curvature of the objective function. H ij = 2 f x i x j = f x i x j Find the numerical values for the matrix H and the vector c corresponding to the cost function f(x), above. (c) Is H always symmetric (H ij = H ji ), (H = H T )? (d) Show that in this example (as in all quadratic problems) the vector is the same as Hx + c. f/ x 1 f/ x 2 f/ x 3 (e) Compute the values of x 1, x 2, x 3 that minimize f([x 1, x 2, x 3 ]). This is an unconstrained minimization. What is the value of the cost function f(x )? You may use Matlab to solve this linear system of three equations and three unknowns. (f) Suppose that in addition to minimizing f, the parameter values must also satisfy two equality constraints: g 1 (x 1, x 2, x 3 ) = 5x 1 1x 2 + 6x = 0 and g 2 (x 1, x 2, x 3 ) = 3x 1 + 9x 2 2x 3 8 = 0 Find the numerical values of the matrix A (2 3) and the vector b (2 1) representing these two equality constraints as Ax b = 0. Show that A ij = g i / x j. (g) Using a column-vector of Lagrange multipliers, λ, and the augmented cost function, f A = f(x) + λ T g(x) = 1 2 xt Hx + c T x + d + λ T (Ax b) (7) CEE 201L. Uncertainty, Design, and Optimization Duke University 2020 H.P. Gavin 7 show that the two conditions f A x = 0 and (x=x, λ=λ ) f A λ = 0 (x=x, λ=λ ) may be expressed in a matrix called the Karush-Kuhn-Tucker (KKT) equation Note that since f A is a scalar H A 0 A T x λ = c b. f A = f T A = 1 2 xt Hx + x T c + d + (x T A T b T )λ (h) Using the numerical values for H, c, A, and b, compute values for the optimal parameter vector x and the Lagrange multipliers λ. You may use Matlab to solve this linear system of five equations and five unknowns. Confirm that Ax b = 0. Are both Lagrange multipliers positive? What is the value of the cost function f(x ) for which Ax b = 0? (i) Change b 2 by 1 and repeat part (h) above. What is the new value of the cost function, f(x )? Compare the difference between f(x ) in part (h) and f(x ) in part (i) to the average of λ 2 in parts (h) and (i). Did you find that f (h) f (i) (λ 2(h) + λ 2(i) )/2? (j) Now without changing b 2, change b 1 by 1 and repeat part (i) above. Did you find that f (i) f (j) (λ 1(i) + λ 1(j) )/2? (k) Are you getting convinced that the values of the Lagrange multipliers indicate the sensitivity of the cost function to changes in the constraints? 8 Homework 3 3. Minimum Total Potential Energy (35 points): The equilibrium configuration of an elastic system minimizes its total potential energy. The total potential energy of an elastic system is the combination of its elastic potential energy (think 1 2 k id 2 i ) and its gravitational potential energy (think mgy j ). In your solids mechanics course you learned how to apply equilibrium equations (think Fx = 0; Fy = 0) to solve for the bar forces in a statically determinate truss. This kind of analysis invokes three assumptions. (a) the truss is statically determinate (b) the truss is stable (c) the displacements of the nodes are small compared to the overall geometry of the truss Maybe you thought that the third assumption shouldn t be warranted. Maybe you thought that the first assumption was too limiting. Maybe you thought that the second assumption couldn t always be tested. Well, the principle of minimum total potential energy does away with all of those assumptions. Pencil-and-paper solutions with this principle are practical for only the most simple problems. However, your nascent mastery of methods of numerical optimization provides the power do away with these kinds of assumptions, allows you to find the displaced configuration of the truss (which the methods of joints does not do), and, frees you from every working through the methods of joints, forever. The five bar truss shown below shows bar numbers, node numbers, initial node locations, and the direction of gravitational acceleration. Is it statically determinate? Don t answer that. By using the principle of minimum total potential energy you don t even need to care. In general each truss bar can have its own stiffness k i = (EA/L) i, Each node has its own coordinates. Coordinates for nodes 1 and 2 are fully pinned: (x 1, y 1 ) = (0, 0), (x 2, y 2 ) = (10, 0). Coordinates (x 3, y 3 ) and (x 4, y 4 ) for nodes 3 and 4 are not. (a) Given values for all coordinate positions, write equations for the length of each bar. (b) Given the initial length of each bar, write equations for the change in length of each bar. (c) Given the stiffness of each bar, write equations for the elastic strain energy in each bar. (d) Given the mass of nodes 3 and 4, a value for gravitational acceleration, and the vertical coordinates y 3 and y 4, write equations for the gravitational energy energy in each node. CEE 201L. Uncertainty, Design, and Optimization Duke University 2020 H.P. Gavin 9 (e) Write a Matlab function five bar analysis.m that has as its first lines... 1 function [f, g] = five_bar_analysis ( param, cnstnt ) 2 % your i n f o here Plots = 1; % 1 : draw p l o t s ; 0 : don t 5 6 k = cnstnt (1:5); % s t i f f n e s s o f each bar, N/m 7 8 Lo = cnstnt (6:10); % i n i t i a l l e n g t h o f each 1, m 9 10 x1 = cnstnt (11); % x c o o r d i n a t e o f node 1, m 11 y1 = cnstnt (12); % y c o o r d i n a t e o f node 1, m 12 x2 = cnstnt (13); % x c o o r d i n a t e o f node 2, m 13 y2 = cnstnt (14); % y c o o r d i n a t e o f node 2, m m = cnstnt (15:16); % mass o f nodes 3 and 4, kg 16 g = cnstnt (17); % g r a v i t a t i o n a l a c c e l e r a t i o n, m/ s ˆ x3 = param (1); % x c o o r d i n a t e o f node 3, m 19 y3 = param (2); % y c o o r d i n a t e o f node 3, m 20 x4 = param (3); % x c o o r d i n a t e o f node 4, m 21 y4 = param (4); % y c o o r d i n a t e o f node 4, m %... w r i t e nine l i n e s or more to compute t h e bar l e n g t h s, t h e bar s t r e t c h e s, 24 % t h e e l a s t i c p o t e n t i a l energy, t h e g r a v i t a t i o n a l p o t e n t i a l energy 25 % and t h e t o t a l p o t e n t i a l energy... and with its last lines g = -1; % t h e r e are no a c t i v e c o n s t r a i n t s in t h i s problem i f Plots figure (1) 46 c l f 47 hold on 48 plot ([ x1 x3], [y1 y3], -b, LineWidth,4); 49 plot ([ x1 x4], [y1 y4], -b, LineWidth,4); 50 plot ([ x2 x3], [y2 y3], -b, LineWidth,4); 51 plot ([ x2 x4], [y2 y4], -b, LineWidth,4); 52 plot ([ x3 x4], [y3 y4], -b, LineWidth,4); 53 hold off 54 axis ([ ] ); %, equal ) ; 55 drawnow; end 10 Homework 3... and use SQPopt to optimize the minimum energy configuration 1 % f i v e b a r o p t.m 2 % your i n f o here k = [ 400 ; 400 ; 300 ; 200 ; 200 ]; % s t i f f n e s s o f each bar, N/m 5 6 x1 = 0; % f i x e d x c o o r d i n a t e o f node 1, m 7 y1 = 0; % f i x e d y c o o r d i n a t e o f node 1, m 8 x2 = 10; % f i x e d x c o o r d i n a t e o f node 2, m 9 y2 = 0; % f i x e d y c o o r d i n a t e o f node 2, m m = [ 2 ; 5 ]; % mass o f nodes 3 and 4, kg 12 g = 9.81; % g r a v i t a t i o n a l a c c e l e r a t i o n, m/ s ˆ x3 = 3; % i n i t i a l x c o o r d i n a t e o f node 3, m 15 y3 = 4; % i n i t i a l y c o o r d i n a t e o f node 3, m 16 x4 = 6; % i n i t i a l x c o o r d i n a t e o f node 4, m 17 y4 = 4; % i n i t i a l y c o o r d i n a t e o f node 4, m % i n i t l a l l e n g t h o f each bar Lo = sqrt ( [ (x3 -x1 )ˆ2 + (y3 -y1 )ˆ2 ; 21 (x4 -x1 )ˆ2 + (y4 -y1 )ˆ2 ; 22 (x3 -x2 )ˆ2 + (y3 -y2 )ˆ2 ; 23 (x4 -x2 )ˆ2 + (y4 -y2 )ˆ2 ; 24 (x4 -x3 )ˆ2 + (y4 -y3 )ˆ2 ] ); cnstnts = [ k ; Lo ; x1 ; y1 ; x2 ; y2 ; m ; g ]; params = [ x3 ; y3 ; x4 ; y4 ]; 29 params_min = [ 0 ; -8 ; 0 ; -8 ]; 30 params_max = [ 12 ; 8 ; 12 ; 8 ]; % a l g o r i t h m i c parameters % d i s p l a y tolx t o l F tolg MaxEvals Penalty Exponent mmax errf 35 options = [ e ]; [ params_opt, f_opt, g_opt, cvg_hst ] = SQPopt ( five_bar_analysis, params, params_min, params_max, options, cnstnts ); plot_cvg_hst ( cvg_hst, params_opt, 100); (f) Minimize the total potential energy for m 3 = 2 kg and m 4 = 5 kg. What are the final node positions of nodes 3 and 4? (g) Repeat for m 3 = 4 kg and m 4 = 10 kg. (double the loads). Did you expect the displacements to double? Did they? Now what are the final node positions of nodes 3 and 4? What s going on here?? To hand in: your written responses to (a) through (d), your printout of your version of your written responses to (f) and (g). five_bar_analysis.m

Related Search

We Need Your Support

Thank you for visiting our website and your interest in our free products and services. We are nonprofit website to share and download documents. To the running of this website, we need your help to support us.

Thanks to everyone for your continued support.

No, Thanks