Description

Classical Optimization Methods

Summary

This laboratory exercise is scheduled for three weeks. The first week is a tutorial meeting that tutors will provide a revision on the topic and guidelines to complete this exercise.

The exercise includes a program based implementation of line search and gradient optimization methods including Newton’s method and the steepest descent method. The implementation details of these methods will be investigated and their characteristics and/or limitations are to be revealed and compared. This exercise also contains the development of a graphical search method to obtain the optimal setting of simple mechanical system parameters. The last part is about the use of the Matlab optimization toolbox to solve the design of the simple mechanical system.

An individual written report is required for this exercise, which carries 8% (of the total course marks), in additional to the marks for the practical work which is 5%. Report submission is to be uploaded to Moodle. Details of the requirements are given in the Assessment section.

 

Part 1 – Line Search and Gradient based optimization

Exhaustive search is a simultaneous search method in which all the experiments are conducted before any judgment is made regarding the location of the optimum point. This method is very inefficient. The dichotomous search method, as well as the Fibonacci and the golden section methods, are sequential search methods in which the result of any experiment influences the location of the subsequent experiment.

In the dichotomous search, two experiments are placed as close as possible at the centre of the interval of uncertainty. Based on the relative values of the objective function at the two points, almost half of the interval of uncertainty is eliminated.

The principle of gradient-based optimization methods is based on finding the root of the gradient of the objective function. This is because by setting the gradient to zero is a condition for a function to be at its minimum or maximum.

When the objective function is continuous, it can be approximated in a Taylor series, which contains terms of the derivatives of the function. The gradient terms are then extracted, set to zero and solved for the variables in the objective function. The Hessian matrix can be formulated from the gradients and Newton’s method produces the optimum solution in an iterative loop.

While calculating the Hessian may be complicated, the steepest descent method makes use of only the gradient information. In particular, the adjustment on the iterative temporary solution is based on assigning the change as the negative of the gradient.

 

Task 1

1. Let the reduction ratio in the golden search method be 𝜌𝜌. Find the specific choice of the golden section ratio in order to reduce the number of objective function evaluations. For a given search range Δ and the accepted tolerance 𝜏𝜏, show how the number of iterations 𝑁𝑁 can be obtained.
2. Given the objective function 𝑓𝑓(𝒙𝒙) of an optimization problem, assume it is continuous where derivatives exist, give a Taylor series expansion about a point 𝒙𝒙_𝑛𝑛 up to the quadratic (2-nd derivative) term.
3. Obtain the derivative of the series with respect to the infinitesimal change Δx = 𝒙𝒙 − 𝒙𝒙_n, and give the expression in terms of the Hessian matrix.
4. Give the iterative equation to obtain the optimum solution 𝒙𝒙^∗ in the form of Newton’s method.
5. Formulate the steepest descent optimization algorithm in an iterative expression including a control coefficient 𝛼𝛼 on the step size.

Put your answers in the written report.

 

Task 2

1. Given an objective function 𝑓𝑓(𝑥𝑥, 𝑦𝑦) = 10(𝑥𝑥 − 2)² + 𝑥𝑥𝑦𝑦 + 10(𝑦𝑦 − 1)², write a Matlab program to show the function landscape. Find the minimum, with a user-specified tolerance of 𝜏𝜏 = 0.01, by using the Golden section method. You can set the search range as 𝑥𝑥 ∈ [−5,5], 𝑦𝑦 ∈ [−5,5].
2. You can illustrate your program output as:

 

3. Verify your results by calculating the gradients, 𝜕𝜕𝑓𝑓/𝜕𝜕𝑥𝑥, 𝜕𝜕𝑓𝑓/𝜕𝜕𝑦𝑦; then solve for the optimum design variables 𝑥𝑥^∗, 𝑦𝑦^∗, based on the first-order necessary condition, ∇𝑓𝑓 = 0. Verify that the solution is a minimum by the second-order necessary condition ∇²> 0, i.e., the Hessian is positive definite. Put your answer in the written report.

 

                 

Task 3

1. Given an objective function 𝑓𝑓(𝑥𝑥, 𝑦𝑦) = (𝑦𝑦 − 𝑥𝑥²)² + (1 −𝑥𝑥)², write a Matlab program to perform the following. (Hint: the min. is at 𝑥𝑥 = 1, 𝑦𝑦 = 1 and 𝑓𝑓(𝑥𝑥^∗, 𝑦𝑦^∗) = 0).
2. Use symbolic commands to calculate the gradient vector 𝑮𝑮 and the Hessian matrix (𝑯𝑯), and its inverse 𝑯𝑯⁻¹.
3. Generate a grid of 𝑥𝑥 ∈ [−5,5], 𝑦𝑦 ∈ [−5,5] in steps of 0.1, draw a surface plot of the function. Determine if the function is convex or not, verify it with the Hessian matrix characteristics.
4. Use the initial guess of the solution 𝒙𝒙₀ = [−1, −1]^𝑇𝑇, construct an iterative loop for Newton’s method and terminates when the change of solution (with respect to the previous iteration) is less than a user-specified tolerance, e.g. 𝜏𝜏 = 0.01.
5. Plot the solutions on the objective function surface and a plot against the iteration count. Your plots may look like the figures below.

 

6. Change the objective function to 𝑓𝑓(𝑥𝑥, 𝑦𝑦) = (𝑥𝑥 − 1)² + (𝑦𝑦 − 2)² and search for the optimum solution again. The resultant plots, from an initial guess at 𝑥𝑥 = −3, 𝑦𝑦 = −3, may be as shown below.

 

7. Observe that the solution is obtained after only one iteration (excluding the initial solution), i.e., iteration 2 as shown). Explain why it is possible. (Hint: inspect the form of the objective function, gradient matrix and/or the Hessian matrix.) Put your answer in the written report.

 

Part 2  Mechanical System Design – 1

The shear stress induced along the z-axis when two spheres are in contact with each other is given by

𝜏𝜏_{𝑧𝑧𝑧𝑧} = 1                   3            ² − (1 + 𝜈𝜈) 1 − 𝑎𝑎𝑧𝑧 tan−1 𝑧𝑧/1𝑎𝑎

𝑝𝑝_{𝑚𝑚𝑚𝑚𝑧𝑧}           2 2 1 + 𝑧𝑧

𝑎𝑎

where 𝑎𝑎 is the radius of the contact area and 𝑝𝑝_{𝑚𝑚𝑚𝑚𝑧𝑧} is the maximum pressure developed at the centre of the contact area, and

1−𝜈𝜈21 1−𝜈𝜈22 1/3

3𝐹𝐹

𝑎𝑎 =                          , 𝑝𝑝𝑚𝑚𝑚𝑚𝑧𝑧 = 𝜋𝜋𝑚𝑚2

𝑑𝑑1 𝑑𝑑2

where 𝐹𝐹 is the contact force, 𝐸𝐸₁ and 𝐸𝐸₂ are Young’s moduli of the two spheres, 𝜈𝜈₁ and 𝜈𝜈₂ are Poisson’s ratios of the two spheres, and 𝑑𝑑₁ and 𝑑𝑑₂ the diameters of the two spheres. In many practical applications, such as ball bearings, when the contact load 𝐹𝐹 is large, a crack originates at the point of maximum shear stress and propagates to the surface, leading to a fatigue failure. To locate the origin of a crack, it is necessary to find the point at which the shear stress attains its maximum value. Formulate the problem of finding the location of maximum shear stress for 𝜈𝜈 = 𝜈𝜈₁ = 𝜈𝜈₂ = 0.3. Then we have

𝑓𝑓(𝜆𝜆) = 0.75 + 0.65𝜆𝜆 tan₋₁ 1 − 0.65

𝜆𝜆

where 𝑓𝑓 = 𝜏𝜏_{𝑧𝑧𝑧𝑧}/𝑝𝑝_{𝑚𝑚𝑚𝑚𝑧𝑧} and 𝜆𝜆 = 𝑧𝑧/𝑎𝑎, and the objective is to maximize 𝑓𝑓(𝜆𝜆).

 

 

 

 

 

 

 

Task 4

Write a Matlab program to implement the following.

1. Use the ‘bracketing’ method to find the search ranges for Golden Section and

Fibonacci methods. (Hing: start with a small 𝜆𝜆, e.g. 0.01).

2. Implement the Fibonacci Section Search method to obtain 𝜆𝜆 where it attains the maximum 𝑓𝑓(𝜆𝜆), use a tolerance of 𝜏𝜏 = 0.01.
3. Implement the Bisection Search method to obtain 𝜆𝜆 where it attains the maximum

𝑓𝑓(𝜆𝜆), use a tolerance of 𝜏𝜏 = 0.01.

4. Your program may produce the following figures.

 

5. Compare the number of iterations needed for these methods, and the maximum values obtained. Comment on their overall implementation complexities. Put your discussion in the written report.

 

Mechanical System Design – 2

Consider the design a uniform column of tubular section, with hinge joints at both ends, to carry a compressive load P=2500 kg_f for minimum cost. The column is made up of a material that has a yield stress (σ_y) of 500 kg_f/cm², modulus of elasticity (E) of 0.85×10⁶ kg_f/cm², and weight density 𝜌𝜌 of 0.0025 kg_f/cm³. The length of the column is 250 cm. The stress induced in the column should be less than the buckling stress as well as the yield stress. The mean diameter of the column is restricted to lie between 2 and 14 cm, and columns with thicknesses outside the range 0.2 to 0.8 cm are not available in the market. The cost of the column includes material and construction costs and can be taken as 5𝑊𝑊 + 2𝑑𝑑, where 𝑊𝑊 is the weight in kilograms force and 𝑑𝑑 is the mean diameter of the column in centimetres.

 

The design variables are the mean diameter 𝑑𝑑 and tube thickness 𝑡𝑡,

𝒙𝒙 = [𝑥𝑥₁, 𝑥𝑥₂]^𝑇𝑇 = [𝑑𝑑, 𝑡𝑡]^𝑇𝑇.

The objective function to be optimized is

𝑓𝑓(𝒙𝒙) = 5𝑊𝑊 + 2𝑑𝑑 = 5𝜌𝜌𝜌𝜌𝜌𝜌𝑑𝑑𝑡𝑡 = 9.82𝑥𝑥₁𝑥𝑥₂ + 2𝑥𝑥₁.

The behaviour constraints can be expressed as stress induced be less than yield stress, and stress induced be less than buckling stress. The induced stress is given by

𝑃𝑃        2500

𝜎𝜎_𝑖𝑖 = =

𝜌𝜌𝑑𝑑𝑡𝑡 𝜌𝜌𝑥𝑥₁𝑥𝑥₂

The buckling stress for a pin-connected column is given by

𝜌𝜌²𝐸𝐸𝐸𝐸 1

𝜎𝜎𝑏𝑏 = 𝑙𝑙2 𝜌𝜌𝑑𝑑𝑡𝑡

where the second moment of area of the cross section of the column is

𝜋𝜋8                 2 + 𝑥𝑥22)

𝐸𝐸 =𝑥𝑥1𝑥𝑥2(𝑥𝑥1

The side constraints are given by

2 ≤ 𝑑𝑑 ≤ 14, 0.2 ≤ 𝑡𝑡 ≤ 0.8

The design objective is to minimize the cost, 𝑓𝑓(𝒙𝒙), such that the constraints in stresses 𝜎𝜎_𝑖𝑖, 𝜎𝜎_𝑏𝑏, and sides 𝑑𝑑, 𝑡𝑡 are satisfied. These include

𝑔𝑔

𝜌𝜌𝑥𝑥1𝑥𝑥2

2500   𝜌𝜌²(0.85 × 10⁶)(𝑥𝑥₁² + 𝑥𝑥₂²)

𝑔𝑔₂(𝒙𝒙) = 𝜌𝜌𝑥𝑥₁𝑥𝑥₂ −              8 × 250₂                     ≤ 0.

 

 

 

 

 

Task 5

Write a Matlab code to carry out an optimum design of the given mechanical system by use of a graphical/search approach.

The search method is based on calculating the objection function over a set of grid points. Then discard those points which do not satisfy the constraints. Finally, the solution is determined by identifying the grid point that gives the smallest objective function value.

1. Define the number of grid points, and specify the ranges according the size constraints.
2. Generate a set of grid points (use the ‘meshgrid’ command), calculate and plot the objective function surface with its contour.
3. Calculate the constraints impose by the limits on the stress.
4. Identify and discard the grid points that do not satisfy the constraints (use the ‘find’ command and set the identified grid point surface to non-numeric, i.e., ‘NaN’).
5. Re-plot the objective function surface on another graph.
6. Find the minimum among the evaluated grid point values (use the ‘min’ command with the minimum value and indices as outputs, and convert to grid points by the ‘ind2sub’ command).
7. Annotate the surface plot with the solution (use the ‘text’ and ‘sprintf’ commands).
8. Evaluate the effect of the number of grid points, hence resolution, on the accuracy of the result.