How to solve this error when I use fmincon?

In the script below, I try to optimize the resiliency of a complex network using fmincon. But after some time of running my script I got the following error:

Operands to the || and && operators must be convertible to logical scalar values.

Error in barrier

Error in fmincon (line 813)

[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

Error in Resiliency_optimization (line 10)

[x,fval,exitflag,output] = fmincon(@objectivefun,x0,A,b,Aeq,beq,lb,ub,@nonlco);

A = [];
b = [];
Aeq = [];
beq = [];
lb =[];
ub = [];
t0 = [ 0.15, 0.2, 300, 0.3];
[x,fval,exitflag,output] = fmincon(@objectivefun,t0,A,b,Aeq,beq,lb,ub,@nonlco);

function f = objectivefun(t)

load('G.mat');
load('Load');
p = t(1);
Trigger = t(2);
bgt = t(3);
alpha = t(4);
N = numnodes(G);
iter = 1;
tmax = 80;
M_iter = zeros(iter, tmax);
for jj = 1:iter
    
    [G_dmg,~,~, needRemoveNode,LoadneedDist,M] = Load_initial(G,8,0,350,0.2,Load);
    
    
    
    alreadyCalled = false; % Say that decision and implement have not been called yet.
    for tt = 3: tmax
        
        
        if any(needRemoveNode)|| ~any(needRemoveNode)
            
            [G_dmg,LoadneedDist,needRemoveNode,M] = Load_Stages(G_dmg,needRemoveNode,LoadneedDist,M);
            
            
            
        end
        if  (M(end) >= (Trigger * N)) && ~alreadyCalled % triggering level
            % Calling for the very first time if we get here.
            [G_dmg, Nodes_p] =  Loads_decision(G_dmg);
            [G_dmg] = Loads_implement(G_dmg, Nodes_p, p, bgt, alpha);
            alreadyCalled = true; % Set flag to indicate we've called it.
        elseif alreadyCalled
            % Once they've been called already, call them regardless of the "if" test.
            [G_dmg, Nodes_p] =  Loads_decision(G_dmg);
            [G_dmg, M] = Loads_implement(G_dmg, Nodes_p, M, p, bgt, alpha);
        end
    end
    M_iter(jj,:)=M;
    %    S = std(M_iter);
    %    SE = S/sqrt(size(M_iter,1));
end

Mavg = mean(M_iter,1);

[~, maxidx] = max(Mavg);
f = find( (Mavg <= 0.1 * N) & (maxidx <= 1:length(Mavg)), 1 );

end

function [c, ceq] = nonlco (t)
c = [];
ceq = [];
end

Thanks

[Matt J]

Your objective function looks quite non-smooth and discontinuous, which violates the assumptions of the fmincon algorithms.

[Waseem AL Aqqad]

Thanks, Matt J. You are right, that is why Walter recommended ga(), particleswarm() or simannealbnd() instead of fmincon.

Answer 1

[Walter Roberson]

        if  (M(end) >= (Trigger * N)) && ~alreadyCalled % triggering level
            % Calling for the very first time if we get here.
            [G_dmg, Nodes_p] =  Loads_decision(G_dmg);
            [G_dmg] = Loads_implement(G_dmg, Nodes_p, p, bgt, alpha);
            alreadyCalled = true; % Set flag to indicate we've called it.
        elseif alreadyCalled
            % Once they've been called already, call them regardless of the "if" test.
            [G_dmg, Nodes_p] =  Loads_decision(G_dmg);
            [G_dmg, M] = Loads_implement(G_dmg, Nodes_p, M, p, bgt, alpha);
        end

What if that condition is never true within a tt value? You would leave M unchanged. Is that appropriate? That is going to affect the mean()

        f = find( (Mavg <= 0.1 * N) & (maxidx <= 1:length(Mavg)), 1 );

What if no entries match that condition? Then f would be empty, and that would cause problems in the optimization function.

You are returning an index. Indices are integers. fmincon() will typically give up easily when it sees integer values, deciding that the function is flat. Your function being minimized should be continuous, not discrete.

[Waseem AL Aqqad]

What if that condition is never true within a tt value? You would leave M unchanged. Is that appropriate? That is going to affect the mean().

Ans: True, M will be unchanged as I'm simulating the delay response of the recovery process (Load_decision and Load_implement functions will not be called for some tt steps).

What if no entries match that condition? Then f would be empty, and that would cause problems in the optimization function.

Ans: You are right. I debugged the script and in some iterations, f was empty. Especially, when the time frame is not wide enough. Any suggestions please?

Any recomendation for designing a continuous function instead of a discrete one. This script is all about the simulation of attacking a graph and recovering it during a time frame tt = 1:tmax. The output M is the number of failed nodes. length(M) is used as time steps. So, If we plot M against the time steps it would look like a plateau( I attack and recover within tt time frame). And eventually I want to optimize the minimum time needed for my recovery algorithm to heal/recover 90% of the graph nodes ( a resiliency metric).

Many thanks.

[Walter Roberson]

For integer outputs, you should use ga() or particleswarm() or simannealbnd(). At the moment I do not remember whether patternsearch() or surrogateopt() can be used for integer outputs.

[Waseem AL Aqqad]

Thanks Walter.

As you mentioned in your comment, " What if no entries match that condition? Then f would be empty, and that would cause problems in the optimization function." I got an error when f is an empty vector as ga() for example accepts only scalar values.

Here, Mavg should be something like this,

Mavg = [0 8 24 99 230 500 1000 2000 3800 4000 3444 2800 1700 666 332 100 50 0] 

in order not to get an empty value of f and match this condition :

f = find( (Mavg <= 0.1 * N) & (maxidx <= 1:length(Mavg)), 1 );

However, in some tt steps, Mavg might be something like this:

Mavg = [0 8 24 99 230 500 1000 2000 3800 4000 4000 4000 4000 4000 4000 4000 4000 4000]

and this will not match the condition. The reason that I get such an Mavg is that I introduce some delay and don't call both decision and implement functions at some tt values:

        if  (M(end) >= (Trigger * N)) && ~alreadyCalled % triggering level ( delay response of recovery process)
            % Calling for the very first time if we get here.
            [G_dmg, Nodes_p] =  Loads_decision(G_dmg);
            [G_dmg] = Loads_implement(G_dmg, Nodes_p, p, bgt, alpha);
            alreadyCalled = true; % Set flag to indicate we've called it.
        elseif alreadyCalled
            % Once they've been called already, call them regardless of the "if" test.
            [G_dmg, Nodes_p] =  Loads_decision(G_dmg);
            [G_dmg, M] = Loads_implement(G_dmg, Nodes_p, M, p, bgt, alpha);
        end

Is there an approach to tackle this issue where I can optimize f when I introduce a delay (Trigger * N).

Thanks again.

[Walter Roberson]

You can use isempty(f) to determine a different value to return insted of empty.

[Waseem AL Aqqad]

Like any random integer? Or should it be chosen wisely?

Something like this?

f = find( (Mavg <= 0.1 * N) & (maxidx <= 1:length(Mavg)), 1 );
if isempty(f)
f = randi(20)
end

Would it affect the optimization process?

[Walter Roberson]

Is it a situation that you want to reject as being unsuitable? If so return a value higher than any valid value.

If it is a situation that indicates that you have found a best-possible configuration, then return 0.

Do not return a random number: that would result in the situation being preferred to some (but not all) valid configurations.

[Waseem AL Aqqad]

Okay then its option 1 for my case.

Thank you so much for your help. I am really grateful.