while loop with user defined function won't cycle through all conditions of the function.

I am utilising a while loop to determine the height of water in a dam and the discharge rate that is dependant on that height. The function that determines the discharge rate has a condition of if the water height is less than 20% then the discharge is 0. The issue is that when the height goes above this threshold the discharge rate determined by the function doesn't reflect that and will populate the rest of the array with 0's.

User defined function

function [Q] = O_R_A(H)
%Operating rule A for varying dam heights
%   Determining the usage rates for the dam for different dam heights
H_max = 30.9; %dam height in m (scaled 1:300)
Q_out = 11923.2; %normal output rate
Q = 0;
    if H >= 1*H_max;
        Q = NaN();
    elseif H >= 0.9*H_max;
        Q = 1.5 * Q_out;
    elseif H < 0.9*H_max & H >= 0.2*H_max;
        Q = Q_out;
    else H < 0.2*H_max;
        Q = 0;
    end   
end

and the while loop that utilises the function is

m = 1;
H_A = 16.2; %initial water height
while m < length(Q_in)
    m = m+1;
    Q_A(m) = O_R_A(H_A);
    delta_a(m) = (((Q_in(m-1) - Q_A(m))/f)*30.9); %change in water level
    H_A(m) = H_A(m-1) + delta_a(m); %new water level
end

when the water level rises back above the 20% threshold I would like for the discharge rates to be calculated as per the conditions in the function.

Answer 1

[John D'Errico]

H_A is a vector. While it starts out as the scalar value 16.2, you are growing the vector in length. This is bad for other reasons, but we will not get into that now. Just recognize that your code will be far better if you preallocate H_A to have the final length in the beginning. Thus start with this:

H_A = zeros(size(Q_in));
H_A(1) = 16.2;

But the failure to preallocate H_A properly is not your problem.

The reason your code fails is because a if/elseif construct does NOT process every element of the vector separately!

Your code looks like this:

   if H >= 1*H_max;
     ,,,

   elseif H >= 0.9*H_max;
     ...

But H is a vector here. An if statement does NOT take a different branch some of the time for each element of the vector! (Also, you don't need to put semi-colons on an if statement test.)

Instead, you need to use indexing.

Q = zeros(size(H));
Q(H >= 1*H_max) = NaN;
Q(H >= 0.9*H_max) = 1.5 * Q_out;
Q((H < 0.9*H_max) & (H >= 0.2*H_max)) = Q_out;
Q(H < 0.2*H_max) = 0;

That last line is superfluous, because Q is already zero there, since it started out as zero in all cases anyway. I added it merely for clarity.

[Rhys Webb]

Am i able to use the code

Q = zeros(size(H));
Q(H >= 1*H_max) = NaN;
Q(H >= 0.9*H_max) = 1.5 * Q_out;
Q((H < 0.9*H_max) & (H >= 0.2*H_max)) = Q_out;
Q(H < 0.2*H_max) = 0;

as part of the function) as this is part of the requirements for this task)?

I have changed my function to reflect the code above and it produces the following error

In an assignment A(:) = B, the number of elements in A and B must be the same.

Error in assignment_code (line 38) Q_A(m) = O_R_A(H_A);

[Guillaume]

You need to think a bit more about what you're doing. You also should learn to debug your code. You would see the problem straight away if you stepped through the code looking at what is actually happening.

You're passing a vector to your function, so John showed you how to deal with a vector input. As a result, the function returns a vector. If you want the function to return a scalar, then pass it a scalar (in which case your original function also works).

[Rhys Webb]

I am currently in debugging mode and I suppose I should have said I don't fully understand the error and what it means. At this point when the error occurs Q_A is showing as a 1x2 matrix as is H_A.

I require a vector as my answer (which will contain 155 points once it is all said and done) and to put it simply I don't understand why this error occurs. I am not sure if it is due to the fact that I am exhausted and sleep deprived or if I am not a very good teacher (having been self taught for the most part).

Any help in the right direction would be appreciated.

[Guillaume]

O_R_A(H_A) is a vector

Q_A(m) is a scalar

You can't stuff a vector inside a scalar. That's what the error is telling you.

The simplest fix is to pass the scalar H_A(m) to your function (original or modified, both work with scalars) instead of the whole vector H_A.

This has nothing to do with the final result which will still be a vector.