MAS252
Dr Sam Rigby – sam.rigby@sheffield.ac.uk MAS252 Coursework
The bisector method for optimisation: ‘Finding things of in- terest’
During the remaining computer lab sessions, we will:
Introduce the bisector method and implement it to solve a fairly simple search problem;
Thus, use the while loop function, rather than the for loop;
Gain an understanding of ballistics problems;
Consider a more general optimisation problem using these tools
Introduction
Imagine that you are designing a roof for a stadium using a piece of software. You pass to that software certain constraints (the load that the roof must bear, the maximum cost, the size of the roof, constraints from the Eurocodes on wind and snow loads etc.) as well as the specification for a range of structural elements that can be incorporated into the design (e.g. the smallest size of strut to incorporate).
You press the button and the software starts to design the building. First, it must reject all solutions that don’t bear the load, are too expensive, etc. However, there are still likely to be a very large number of appropriate solutions. You may wish for the software to present you with the “best” one, or to present you with a range of choices of “good” results.
How does the software make these choices? This is the optimisation problem.
If we can form an appropriate metric (e.g. the cost) to evaluate all the results as a function of the input parameters then within the model we may end up with a very complex cost surface (Fig. 1). We may then seek the global minimum on this surface (the “best” result) or select local minima from different parts of the parameter space (the “good” results). However, we need to find these and searching a complex, high dimensional surface would be very computationally expensive. We need a computational method that will quickly “zoom” in to the solution.
The while loop (recap)
So far on the course, we have mainly used the for loop structure. e.g.:
for loop1=1:n
…
end
Note that an equivalent statement would be:
count=1;
while count <= n ... count=count+1; end
1
Figure 1: A hypothetical cost surface as a function of two hypothetical control variables
Clearly, this seems pretty pointless as we are using extra lines of code to do the same thing. However, the while loop is exceptionally powerful when we don’t know how many times we need to repeat our loop.
Consider the problem described in the introduction – we need to keep looping until we find the minimum but we don’t know how long this will take. Thus, an example while loop might say:
while error > tolerance …
end
The Bisector method
This is a simple method for finding a value of interest (such as the minimum of a function). The algorithm proceeds as follows:
1. We first define constraints to the search (the smallest and largest values that we think could contain the answer). We denote these as left and right, respectively.
2. We then define a target position as mid = (left + right)/2 and evaluate the function at all the three positions, to give yL, yM , and yR. If one of these values falls within our error tolerance (i.e. for root finding this would be close to zero) then we have found the solution and we terminate the search. Otherwise, the next step then depends on the values found:
If yR is closest to the solution, right stays the same and mid becomes the new left. We then repeat step 2;
If yL is closest to the solution, left stays the same and mid becomes the new right. Go to step 2;
If yM is closest to the solution then there are various ways to proceed. The simplest is to see which of yR and yL is closest to the solution and use this to guide the choice between the two options already presented.
However, if the true solution is just to the left of mid but right is closer than left we will never find the true solution. Hence, perhaps a better approach is to not move whichever of
2
left and right is closest and move the other in, but rather cautiously, e.g. if mid is closest, followed by left, we adjust right using:
right = mid + 0.8(right-mid)
3. Once the search is terminated, the value returned is reported back to the user.
The task
We will consider the trajectory of a baseball, where all legitimate solutions to the problem are those that reach (or exceed) a target distance. Specifically, our target distance is the size of a baseball field and we are trying to hit a home run. Our boundary is 90 m away, and our projectile may be launched at any angle, θ, between 0 and 90◦.
As well as being able to control the launch angle, we can also control the aerodynamic drag of the ball. A sphere normally has a drag coefficient, Cd = 0.5. However, by imparting a particular spin to the ball, we can reduce this to as low as Cd = 0.4. In addition, we can control the initial velocity imparted to the projectile between the limits of V0 = 45 m/s and V0 = 60 m/s. Our projectile has a mass, m, of 0.145 kg and a cross-sectional area, A, of 42 cm2. The density of air, ρ, is 1.225 kg/m3.
The initial velocity V0 will act in the direction of the launch angle θ. For what follows, this needs to be decomposed into horizontal, x, and vertical, z, components: U0 and W0, respectively.
We will assume that the drag exerted on the projectile by the air is linearly proportional to velocity, which simplifies the problem. Integration of Newton’s second law gives the horizontal velocity as:
U(t) = U e− k t (1) 0m
and the distance travelled (sx = 0 for t = 0):
s (t) = mU (1 − e− k t) (2)
xk0m
where k is the weight of the ball with respect to its terminal velocity (k = mg/vt). The terminal velocity
is a function of the weight, cross-sectional area and drag coefficient:
vt = 2mg (3)
ρCdA
Here, mass (m) is in units of kg, area (A) is in units of m2, and velocities are in units of m/s.
The equations for the vertical motion are a bit more complicated:
kk
s(t)=−mgt+m(W +mg)(1−e−kt) (4)
W(t)=−mg+(W +mg)e−kt 0m
z0m kkk
For simplicity we can assume that the ball is launched from a vertical height of 0 m.
Step 1 (you should aim to complete a step each week)
Write a function called ballistic.m that solves this set of distance equations and returns to the user the horizontal distance travelled by the ball, assuming that the ball is launched at ground level and comes to rest upon impact with the ground.
Hint: How can we tell when the ball has travelled its maximum distance? What information is MatLab storing that will reveal when the ball hits the ground?
The first line for this function should be something like: 3
function distance = ballistic(theta,V0,Cd,mass,area,tstep,plotting)
Hence, the user inputs the angle, initial velocity, drag coefficient, mass and area of the ball, as well as the time step used to solve the distance equations. If plotting has a value of 1, your code should also return a plot of sz against sx using the plot function. Your code should return nothing if plotting̸= 1. Ensure that you select a value of tstep sufficiently small such that results are precise to ±0.05 m. You might want to play around with this value to get a feel for how the accuracy of your solution is influenced by the time-step size.
Step 2
Write a function called bisector maxdist.m that calculates the launch angle that will give the maximum horizontal distance for any given values of Cd and V0.
The function should take input arguments called start and finish, as well as Cd, V0, mass, area and tstep which will be used to feed into our ballistic.m function. At this stage we won’t need a trajectory plot, so when we call up ballistic.m we should set plotting= 0.
The variables start and finish are our first guesses of left and right for use in the bisector method, which will be used in Step 2 to calculate the angle which gives the maximum distance travelled for a given set of starting conditions.
Our error function used to terminate the while loop should be based on the difference in the distance travelled at mid and the mean of left and right, i.e. |yM − (yL + yR) /2| < 0.05 m.
Q: What is the maximum distance travelled for our least favourable conditions (Cd = 0.5, V0 = 45 m/s)?
Step 3
For a given Cd and V0 combination, there will be two points where the θ vs horizontal distance travelled curve crosses 90 m, i.e. there are two values of θ that will give exactly 90 m, with any angle between giving > 90 m travel distance.
We know from Step 2 that bisector maxdist.m will always return a launch angle that gives us a distance travelled greater than 90 m. We therefore know that there will always be one angle greater than this that will achieve an exact distance of 90 m, and one angle less than this that will achieve an exact distance of 90 m.
Write a new function called bisector bestangle.m, adapted from the function written in Step 2, that now calculates the exact angle that will achieve a home run (= 90 ± 0.05 m travel distance). You need to decide which values of start and finish you feed into your bisector.m function to ensure that the upper and lower values of θ (that will give you a home run) are returned. Or, more specifically, you need to determine a way to ensure that MatLab does this for you.
Our error function used to terminate the while loop should be based on the difference in the distance travelled at mid and 90 m, i.e. |yM − 90| < 0.05 m.
Step 4
Write a MatLab code such that the maximum and minimum launch angles (that result in a home run) are calculated for all values of Cd and V0.
You will need to find suitable intervals for Cd and V0, remembering that as a two-parameter problem, you will have to find 2 × n × m solutions, where n is the number of values you select for Cd, m is the number of values you select for V0, and the 2 comes from the fact that both an upper and lower angle are calculated.
The best way to do this is to have for loop structures that loop over changing values of Cd and V0, and feed the appropriate values of start and finish (i.e. the optimal value determined from bisector maxdist.m, selected as either start or finish depending on whether you want the largest or smallest angle) into your bisector bestangle.m function, which is run separately to determine the upper and lower launch angles for each Cd and V0 combination.
It may take your script several minutes to calculate this. It is a good idea to use the disp([i j]) function to display your loop counters (i and j) during the analysis so you know how far it’s got!
4
Step 5
You are the baseball team coach. You know that in order to invest in a player to hit the ball at a given speed (between the 45 m s−1 and 60 m s−1 limits) costs:
Cspeed = 25000 + (17(45 − V0))2 (5) dollars, and that the cost to train them to apply the appropriate spin to change the drag coefficient is: Cspin = 5000(20(0.5 − Cd))4 (6)
dollars.
The greater the range of angles that achieve a home run, the more likely a home run is to be scored. Hence, you should apply a suitable weighting function to the cost, based on to the range of angles (∆θ = θmax − θmin) that give a horizontal distance that exceeds 90 m (this is effectively risk management)
0.1
Cangle = 1 − (0.9(∆θ − 10)/30)
1
when ∆θ > 40◦
when 10◦ ≤ ∆θ ≤ 40◦ (7) when ∆θ<10◦
Q: Calculate the cheapest way to produce an effective home run scoring baseball player, i.e. determine the minimum cost of any Cd and V0 combination, where the total cost is given as (Cspeed + Cspin) × Cangle
5
The coursework
Once you have completed steps 1–5, you can begin the assessed part of this year’s coursework.
Your task is to write a short (4–6 page) commentary on your MatLab code, detailing your understanding of the code and explaining how and why it operates as it does. You will need to explain why you have used certain MatLab commands, and why you have structured your script and functions in the way you have.
This should include a description of your ‘upper level’ MatLab script (the one that calls upon your functions) as well as the functions themselves. Feel free to use diagrams/flow charts to explain your code.
The coursework is purposefully open ended to give you an opportunity to explain your understanding of the code. Whilst I am happy for you to write your code in groups/pairs, this is an individual submission and any unfair means will be treated in accordance with the Universities policy.
As a guide, you might want to discuss:
How the code works
Why you used certain commands, e.g. while loops
Any issues you faced, and how you successfully debugged them
A discussion of the engineering implications of your code, and what you have learned from this exercise (i.e. the need to balance the accuracy of the solution with computational time)
Your lowest-cost values of V0 and Cd, and how you arrived at those values. As an engineer, is it important to have a single, optimal value, or a range of permissible values, given the accuracy to which things can be built?
The marking scheme is as follows:
Your ability to convince the reader of your understanding of the code: 50%
Written quality of the report: 20%
Accuracy of lowest-cost values of V0 and Cd: 10%
Use of diagrams/flow charts to explain concepts: 10%
Understanding of the engineering context of the bisector method/minimisation problem: 10%
The coursework forms the entirety of your assessment for the programming aspect of MAS252, and hence accounts for 20% of your total module mark.
The deadline is 9:00 am on Friday of Week 11.
6