Aqua Phoenix
     >>  Lectures >>  Matlab 7  


7.1 Exercise 3

For this exercise, we will simulate a terrain of some varying elevation, on which a rider drives the Segway until the Energy in the batteries is exhausted. We will then be able to predict how long the batteries will last for this terrain.

The simulated terrain resembles a 2-dimensional curve, which is divided into regular intervals. In each interval, the curve is either has a positive slope, a negative slope, or a slope of zero for flat terrain.

Figure 7.1
Click image to enlarge, or click here to open
In order to figure out at which point the batteries are exhausted, we will need to compute the energy required for each interval. We can then gradually subtract these intervals from the total energy in the batteries until we hit 0 or fall below.

We create a new m-file for this exercise, named exercise_terrain.m. Since we will use most of the variables from exercise_batterypower.m, we include the m-file on the first line.

Figure 7.2
Click image to enlarge, or click here to open
We need to define 2 vectors of equal length, one of which resembles the regular intervals along the x-axis, and the other resembling the height, whether relative or absolute, of the bounds of each interval. We will space the regular intervals on the x-axis 1 km apart, even though we could certainly generate a much finer sampling:

x = [0:1000:20000] The vector of heights will either need to be defined by a random function generator, or by enumerating each point individually:

y = [50, 70, 150, 130, 80, 80, 80, 10, 50, 50, 10, 140, 200, 200, 250, 305, 160, 150, 150, 150, 100] The units for both x and y vectors is set at meters. While it does not affect the solution how elevated the overall terrain is (mountains or sea level), the relative distances between any two consecutive points is important. For example, over a range of 1000 meters, we would rarely observe a 1000 meter height difference, as this would produce a 45° angle (quite steep).

Figure 7.3
Click image to enlarge, or click here to open
It is not a bad idea to plot the terrain before proceeding. This particular terrain is pictured in Figure 7.1.

Usually, the scale of x- and y-axes is adjusted automatically and independently by Matlab. If the range of values on the x-axis is between 0..20000 and the range on the y-axis between 0..500, the resulting graph is obviously not proportional. In the case of our terrain, we would like to see a proportional or close to proportional scale, hence we include the following statement after the plot function:

axis([0 20000 0 2000]) This will set a range of [0..20000] on the x-axis and a range of [0..2000] on the y-axis, thus leaving a scale of 10:1 between x and y. We must note this scale when interpreting the resulting graph.

Figure 7.4
Click image to enlarge, or click here to open
Next, we compute a few other helpful values, such as the slope of each interval (angle θ) and the actual distance ds, which depends on the slope and is thus not constant like the intervals along x. These values are important for the component of Forceslope.

Figure 7.5
We first assign the length of either the x or y vector to a simpler variable n, which now denotes the last element of the vector:

n = length(y) The slope of each interval is computed by dy/dx, which in terms of a given interval means:

[y(i + 1) - y(i)] / [x(i + 1) - x(i)] To do this quickly for the entire vectors of x and y values, we simply take sub-vectors and perform the array operations on them. Specifically, in the term [y(i + 1) - y(i)], y(i + 1) is always one element after y(i). We can thus take a sub-vector of y with values ranging from [2..n] and subtract a second vector with values ranging from [1..n-1]. We do the same with vector x, and finally perform division on the resulting vectors:

s = (y(2:n) - y(1:n-1)) ./ (x(2:n) - x(1:n-1))

Figure 7.6
Click image to enlarge, or click here to open
We can now convert the vector of slopes to angle values in degrees, which later becomes a parameter to Fslope:

alpha = atan(s) * 180/pi

Figure 7.7
Click image to enlarge, or click here to open
We can perform a similar operation on computing the actual distance values ds within each interval (i, i + 1):

ds = √(((y(2:n) - y(1:n-1)) .^ 2) + ((x(2:n) - x(1:n-1)) .^ 2)) which computes for each interval the common equation s = √(x2 + y2).

Figure 7.8
Click image to enlarge, or click here to open
We now proceed to computing the required Energy for each interval. We first re-define variable n to reflect the length of the previously computed vectors, which, due to differentiation, is one less than the original length:

n = length(ds) We also initialize an empty results vector dE:

dE = [ ]

Figure 7.9
Click image to enlarge, or click here to open
We now loop through each element of vector alpha, which enumerates each interval's slope in degrees:

for i=1:n


Figure 7.10
Click image to enlarge, or click here to open
Within the loop we compute the energy used for each interval using equations , , and :

Torque   =   (dwheel / 2) * Ftotal(FSA(hrider, mrider) + Asegway, vkmh2msec(vconst), cD, msegway + mrider, alpha(i))
Pmech effective   =   Torque * v_msd2rad(vkmh2msec(vconst), dwheel)
dE(i)   =   Pmech effective * ds(i) / vkmh2msec(vconst) * (1 / effoverall)
Figure 7.11
Click image to enlarge, or click here to open
We can now evaluate the results vector dE, which enumerates the energy requirements for each interval.

Figure 7.12
Click image to enlarge, or click here to open
We proceed to computing the cumulative energy level, given fully charged batteries with Ebatteries total. For each interval in dE, we will deduct dE(i) from Ebatteries total. We begin by initializing a results vector Ecsum, which will hold the cumulative energy level after each interval:

Ecsum = [ ] Since at the very beginning the batteries are still fully charged, the cumulative sum of energy in the batteries at index 1 equals Ebatteries total:

Ecsum(1) = Ebatteries total

Figure 7.13
Click image to enlarge, or click here to open
We now build a loop that iterates over each interval:

for i=1:n

where the body of the loop will handle the decay of energy.

Figure 7.14
Click image to enlarge, or click here to open
There are two cases we need to handle, and we will simplify the problem to some extent. When driving on even or up inclined terrain (alpha >= 0), we obviously require energy. However, when riding downhill (alpha < 0), we will assume that no energy is required, even though there may be aerodynamic drag. Thus, we merely need to differentiate between positive and negative values of alpha as follows:

if (alpha(i) >= 0)
where statement1 contains the code for required energy, and statement2 contains the code for no required energy.

Figure 7.15
Click image to enlarge, or click here to open
For positive and zero alpha values, we need to deduct the required energy from the previous cumulative sum as follows:

Ecsum(i + 1) = Ecsum(i) - dE(i) For negative alpha values, we simply retain the previous cumulative sum, and we do not add or subtract any energy from the battery:

Ecsum(i + 1) = Ecsum(i) The final loop becomes:

for i=1:n
  if (alpha(i) >= 0)
    Ecsum(i + 1) = Ecsum(i) - dE(i)
    Ecsum(i + 1) = Ecsum(i)


Figure 7.16
Click image to enlarge, or click here to open
We can now evaluate the results vector Ecsum.

Figure 7.17
Click image to enlarge, or click here to open
Finally, we plot the decay of energy in the batteries:

plot(x, Ecsum), xlabel('Distance in meters'), ylabel('Energy left in batteries'), title('Energy used on Simulated Terrain'), grid

Figure 7.18
Click image to enlarge, or click here to open
The resulting graph shows us a decreasing and at times constant amount of Energy in the batteries.

Figure 7.19
Click image to enlarge, or click here to open
While we can gather from the graph approximately when the batteries are empty, we would like to have an exact figure printed in plain English to the screen.

For this endeavour, we will need to iterate over the cumulative sum vector and find the point at which the sum drops below 0. Once we have found this spot, we will need to interpolate between the positive and negative points on the graph to get an exact figure. We begin by creating a new loop. This time, we will use a while loop, which is just another way of iterating over a range.

For a while loop, we initialize a loop variable to some quantity. In the while statement, we check whether this quantity fulfills a certain criteria, and as long as it does, the loop keeps on iterating. In our case, the loop variable will be the index of Ecsum. We initialize it to 2, because we will start by checking the 2nd element of the cumulative sum, knowing that the first element must be positive:

i = 2 We also initialize a variable for the total distance travelled to 0. For every interval, we will then add the amount travelled over that interval:

dtravelled = 0

Figure 7.20
Click image to enlarge, or click here to open
We now write the loop statement:

while (i <= n)
  dtravelled = dtravelled + ds(i - 1)
  i = i + 1
where the code in the loop is iterated over as long as the index into Ecsum is valid, i.e. in the range until and including the last element. The statement dtravelled = dtravelled + ds(i - 1) adds the distance travelled over the past interval. The statement i = i + 1 increments the loop variable. If we did not have this statement, the loop could possible iterate forever and never stop.

Figure 7.21
Click image to enlarge, or click here to open
For each value of Ecsum that we iterate over, we would like to test whether or not it falls below 0. When it does, we would like to aggregate the distance travelled, display it onscreen, and stop the loop.

First, we test for the case in which Ecsum is negative:

if (Ecsum(i) < 0)

where the code before the term break aggregates and displays the distance. The keywork break is used to stop a loop.

Figure 7.22
Click image to enlarge, or click here to open
When we find a negative value of Ecsum, we need to interpolate between the last positive value and the current negative value, in order to find the crossing point of the curve with the x-axis.

For cumulative energy, the current value Ecsum(i) is the first negative value, and the last positive value must be Ecsum(i - 1).

For distance travelled, the current value is x(i), and the last value is x(i - 1).

We now need to interpolate in the range of [x(i) .. x(i - 1)] to find the crossing point of the line with points (x1, y1) = (x(i - 1), Ecsum(i - 1)) and (x2, y2) = (x(i), Ecsum(i)). Matlab has a built-in function for interpolation, named interp1. Use the help function to find out specifics on how this function works: help interp1.

xcross = interp1([Ecsum(i - 1), Ecsum(i)], [x(i - 1), x(i)], 0)

Figure 7.23
Click image to enlarge, or click here to open
Once we have the point of intersection with the x-axis, we can compute the fraction of distance in the current interval, for which the energy level in the batteries is still non-zero:

dremain = ds(i) * ((xcross - x(i - 1)) / (x(i) - x(i - 1))) Finally, we add this part of the interval to the total distance travelled, and convert to km by dividing by 1000:

dtravelled = (dtravelled + dremain) / 1000

Figure 7.24
Click image to enlarge, or click here to open
To print the results in plain English, we use a function from the C programming language called sprintf:

disp(sprintf('The batteries will be empty after %f km.', dtravelled)) where %f is a placeholder for floating point numbers. (%d is the placeholder for integers)

Figure 7.25
Click image to enlarge, or click here to open
We finally evaluate exercise_terrain.m in the Command Window and obtain the graph of decreasing cumulative energy in the batteries, and a verbal answer as to when the Segway's batteries are drained.

Figure 7.26
Click image to enlarge, or click here to open