Finding maximum of a function using optimization toolbox

function [p_uth,Eload] = threshold(ann_load)
for i = 1:size(ann_load,2)
    p_uth(i) = 1.6; % upper threshold - temporarily assigned for now
    pES_val = ann_load(:,i)-p_uth(i); % calculating load values above the upper thresold value
    load_val = ann_load(:,i);  % assigning load values of a day to load_val variable
    x_p_uth = findX((1:size(ann_load,1))',load_val,p_uth(i)); % finding x-coordinates where upper threshold touches the load curve 
    y_temp = zeros(size(x_p_uth,1),1)*p_uth(i); % array of upper thresold values - to facilitate calculations
    pES = pES_val(pES_val>0);  % P_ES values
    index_pES = find(pES_val>0);  % index of P_ES values greater than 0
    y_val = [y_temp;pES]; % uniting arrays of y values
    x_val = [x_p_uth;index_pES]; % uniting arrays of x values
    [x_sorted,sort_index] = sort(x_val); % sorting the x values array
    y_sorted = y_val(sort_index); % sorting the y values array
    Eload(i) = trapz(x_sorted,y_sorted)*5/12; % finding area under the curve of pES values
end
end

In the above I need to find the maximum of 'Eload', which is the objective function. As evident from the code, I am considering p_uth and pES as the variables. Hence Eload = f(p_uth,pES). I need to maximize Eload which can be done by varying p_uth, which needs to be varied between its upper and lower bounds. Eload and p_uth will be 1x366 array. The input to the function, ann_load, is a 288x366 array.

pES is calculated automatically based on the relation mentioned in the above code

pES_val = ann_load(:,i)-p_uth(i); 
pES = pES_val(pES_val>0); 

It can be noted that array size of pES varies in every iteration as the no. of positive values change in every iteration, based on the above equation.

Given the current conditions, I don't know how to use the optimization toolbox to find the maximum. I think I can solve it using nested loops in a scripted file instead of a function but I want the best possible solution using optimization toolbox. Does anybody have an insight as to how to solve it. Any help is appreciated.

[Matt J]

In the above I need to find the maximum of 'Eload', which is the objective function....Eload and p_uth will be 1x366 array.

If Eload is vector-valued, would does it mean to "maximize" it? Are you solving a separate optimization problem for each i?

[Satya Venkata Siddhardha Manchala]

Yes. ann_load is a 288x366 array i.e. each of the 366 columns give a day's data. The data is provided in 5-minute intervals hence there are 288 readings. Eload should be a 1x366 array because Eload has the maximum value of each day for all 366 days.

Answer 1

[Matt J]

Assuming you are solving a separate optimization for each i (see my comment above), it seems to me that the solution is very simply

p_uth=min(ann_load)

This will maximize the elements of pES_val and hence its integral.

[Satya Venkata Siddhardha Manchala]

Sorry for not providing this information earlier. The variables pES, p_uth and the objective function Eload have bounds.

 0 < pES <= max_cap 
 0 < p_uth < max(ann_load) 
 0 < Eload <= tot_cap

Since pES has bounds between 0 and max_cap, p_uth cannot be min(ann_load).

I hope this makes sense.

[Matt J]

Regardless, the objective Eload is monotonic in p_uth. The smallest p_uth value that satisfies the bounds will therefore have to be the optimum.

[Satya Venkata Siddhardha Manchala]

The value of p_uth changes everyday since ann_load(:,i) changes with i. Hence every value in the 1x366 array is unique. As you have said, the idea is to find the smallest value of p_uth such that pES and Eload are maximum while staying within their respective bounds. I tried to do that by applying fmincon using the syntax

x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)

But as you can see, I have to give the starting point x0 of both variables. I could give starting value for p_uth as min(ann_load) as its size is 1x366 through out the iterations whereas I cannot say the same for pES since its size varies in every column. Similary I needed to give the bounds of the objective Eload somehow to the toolbox. So I created the constraints function

function c = eqcon(p_uth)
c = Eload-13.5;
end

Although I was pretty sure that I was wrong, I still needed to input the constraint on Eload.

Now since I did not know how to specify pES to fmincon, I ignored pES and considered only p_uth as the variable. Now I used the above syntax as

[x,fval]=fmincon(@threshold,min(ann_load),[],[],[],[],zeros(1,366),max(ann_load),@eqcon)

When I executed the above command, Matlab gave the error "Supplied objective function must return a scalar value." Now I don't know how to proceed further. Can you help me get out of this ?

[Matt J]

Well, you cannot use fmincon, because fmincon assumes that your objective function is differentiable as a function of your unknowns which is not the case here.

However, I do not know why you are pursuing fmincon or any numerical solver for that matter when, as I explained to you in my last comment, we already know the solution, even when bounds on p_uth are present:

p_uth(i)=max( min(ann_load(:,i) ),0 )

[Satya Venkata Siddhardha Manchala]

Actually p_uth would not always be close to min(ann_load(:,i)). Sometimes it would be close to max(ann_load(:,i)) as well. This depends upon the difference between maximum and minimum of ann_load values of the day. In cases with a much bigger difference, your formula finds p_uth value within the bounds, but it will produce a pES value that is way out of its bounds. This can be better explained with data.

Looking at the data that I have, max(ann_load) - min(ann_load) value ranges from 0.154 to 17.955. Let me consider the extreme values for demonstration.

0.154 occurs for i = 97 and 17.955 occurs for i = 236.

For i =97,

max(ann_load(:,97)) = 0.4134
min(ann_load(:,97)) = 0.2594

Using the formula you mentioned, I would get p_uth(97) = 0.2594. This would give me a 287x1 array as output for pES because there are 287 values above the minimum. Since max(pES) = 0.154 < max_cap which is 5 and Eload = 9.5141 < tot_cap which is 13.5, output is within bounds and your formula worked for i = 97.

Now let us consider for i = 236

max(ann_load(:,236)) = 17.9901
min(ann_load(:,236)) = 0.0351

Using your formula, I'd get p_uth(236) = 0.0351. Similar to results for i = 97, I'd get a 287x1 array for pES. But in this case, max(pES) = 17.955 > 5 and Eload = 317.4562 >> 13.5. So in this case, p_uth should not go below 12.9901 as that would make pES go beyond its upper limit.

Hence depending upon the day i.e. i value, p_uth should change. Also, I don't think I have mentioned this earlier but p_uth value doesn't necessarily have to be a value from the 288 readings in each column. It just has to be a value between minimum and maximum readings in the column. Hope this helps you in better understanding my predicament.

[Matt J]

Okay. Well, we may have to do this numerically, but fminbnd is better for this than fmincon. fmincon is overkill for a problem with a single unknown, and has technical differentiability assumptions that are violated here, as I mentioned.

function Main
    
  max_cap=...
  tot_cap=...    

    for i=1:size(ann_load,2)
        
        p_uth(i)=isolve(ann_load(:,i));
        
    end
    
    function p_uth_Optimal=isolve(ann_load)
    
    %% Derive bounds
    
    lb=max(max(ann_load-maxcap),0); %lower bound on p_uth
    ub=max(ann_load); %upper bound on p_uth
    
    p_uth_Optimal=fminbnd(@objective,lb,ub);
    
        function Eload=objective(p_uth)
            
            Eload= trapz(max(ann_load-p_uth,0)); 
            
            if Eload>tot_cap, 
                Eload=inf;
            else
                Eload=-Eload; %convert maximization to a minimization
            end
            
        end
    
    
    end

end

[Satya Venkata Siddhardha Manchala]

That worked perfectly! Thank you for doing it.

Yes, I understand now that fmincon is not at all suitable for my problem. At the time I was formulating the problem, I did not know how to add constraints on the objective Eload. So I thought adding it as an inequality constraint in fmincon would solve it. It did not strike me that fminbnd could be used to find the optimal p_uth value.

Can you explain what you were doing in the 'if' loop after finding Eload value using trapz. I assume that Eload = -Eload because fminbnd finds minimum value and we need the maximum, which can be done by putting a -ve sign infront of Eload. But what about Eload = inf. What would that achieve ?

[Matt J]

That worked perfectly! Thank you for doing it.

You're welcome, but please Accept-click the answer to indicate that the problem is resolved.

But what about Eload = inf. What would that achieve ?

That is my way of enforcing Eload<tot_cap. By setting the objective to inf in the region that violates this constraint, we communicate to fminbnd that the minimum is not located there.

[Satya Venkata Siddhardha Manchala]

That's genius! That thought would have never occured for me. I guess having a ton of experience like you will help in solving problems with the least amount of effort. Any tips on how to get better at writing code?

[Matt J]

Any tips on how to get better at writing code?

Stay active in the Answers forum... maybe also view Mathworks webinars from time to time...

https://www.mathworks.com/videos.html

[Satya Venkata Siddhardha Manchala]

Thanks a lot for the info and for solving my problem. I'll take your advice and try to improve my coding skills.