Genetic algorithm optimization tool not accessing the objective function

I am trying to optimise an objective function with the name objective function. Optimisation parameter is Charge/Area and Frequency is meant as an constraint which is in objective function itself with an if condition.

I am using Parallel computing. But I observed parallel pool is left idle. Then, I put a breakpoint at the first line of objective function but the code never paused indicating there is something wrong with the code. I would be glad if somebody can suggest what went wrong in this.

clc;
clear; 
close all;
%% Initial Conditions
%starting point
L_b_0=200;
w_b_0=60;
L_m_0=200;
w_m_0=400;
t_hSi_0=98;
t_aln_0 = 477;
L_s_0 = 150;
gap_s = 5;
w_s_0 = 0.5*(w_m_0 - w_b_0 - 2*gap_s);
x0=[L_b_0 w_b_0 L_m_0 w_m_0 t_hSi_0 t_aln_0 L_s_0 gap_s]

x0 = 1×8

200 60 200 400 98 477 150 5

objective(x0)

ans = 100000

%% Bounds and Constraints
%bounds
lb=[100 5 100 5 10 100 100 5];
ub=[500 150 500 1000 350 1000 1000 800];
%linear inequalities
A=[-0.3 1 0 0 0 0 0 0; 0 1 0 -1 0 0 0 2];
b=[0; -5];
%linear equalities equalities
Aeq=[];
beq=[];
%nvars
nvars=8;
% CalcCharge_RB([200,156.3643632,200,360.5678223,239.2755624])   
% Constraint_RB([200,156.3643632,200,360.5678223,239.2755624])
%Call solver to minimize objective function
%% Optimisation
opts = optimoptions('ga','OutputFcn',@Anju_ga_save_each_gen_NonLin,'InitialPopulationMatrix',x0,'PopulationSize',5000,'EliteCount',20,'CrossoverFraction',0.4,'HybridFcn','fmincon','UseParallel',true);
xopt = ga(@objective,nvars,A,b,Aeq,beq,lb,ub)%,[],opts)

ga stopped because the average change in the fitness value is less than options.FunctionTolerance.

xopt = 1×8

500 5 100 1000 10 100 100 5

objective(xopt)

ans = 100000

% Charge = CalcCharge_RB(xopt);
%Freq = Constraint_RB(xopt);
%% Objective Function
%objective function - Charge/Area + Constraint - Frequency
function obj = objective(x)
%Defining Parameters
L_b=x(1)*10^-6;
w_b=x(2)*10^-6;
L_m=x(3)*10^-6;
w_m=x(4)*10^-6;
t_hSi=x(5)*10^-6;
t_aln = x(6)*10^-9;
L_s = x(7)*10^-6;
gap_s = x(8)*10^-6;
w_s = 0.5*(w_m - w_b - 2*gap_s);
t_si = 1*10^-6;
t_sio2 = 3*10^-6;
t_al=0.1*10^-6;
if(w_b <= 0.3.*L_b && w_s >= 5)
%% Calculating Required parameters
%Densities
rho_al = 2700;
rho_aln = 3300;
rho_si = 2329;
rho_sio2 = 2200;
%Elastic Modulus of Materials
E_al = 70*10^9;
E_aln = 330*10^9;
E_si = 188*10^9;
%neutral axis
n = (E_al.*t_al.^2 + 2.*E_al.*t_al.*t_aln + 2.*E_al.*t_al.*t_si + E_aln.*t_aln.^2 + 2.*E_aln.*t_aln.*t_si +E_si.*t_si.^2)/(2.*(E_al.*t_al + E_aln.*t_aln + E_si.*t_si));
I_al = (1/12).*w_b.*t_al.^3 + w_b.*t_al.*(t_aln + t_si - n + 0.5.*t_al)^2;
I_aln = (1/12).*w_b.*t_aln.^3 + w_b.*t_aln.*(t_si - n + 0.5.*t_aln)^2;
I_si = (1/12).*w_b.*t_si.^3 + w_b.*t_si.*(0.5.*t_si - n)^2;
m = L_m.*w_m.*(rho_aln.*t_aln + rho_si.*t_si + rho_sio2.*t_sio2 + rho_si.*t_hSi);
Force = m.*9.8;
d_31 = -1.9159*10^-12;
m_s = L_s.*w_s.*(rho_aln.*t_aln + rho_si.*t_si + rho_sio2.*t_sio2 + rho_si.*t_hSi);
%% Charge Calculation
Charge = (d_31.*9.8.*w_b.*L_b.*t_aln./I_aln).*((0.5.*m+m_s)*L_b - m_s.*L_s + 0.5.*m.*L_m);
%% Area Calculation
Area = (L_b+L_m).*w_m;
%% Figure of Merit
Normalised_Charge = Charge.*10^13;
Normalised_Area = Area.*10^9;
Fom = Normalised_Charge./Normalised_Area;
obj = Fom;
%% Frequency
E_tot= E_al.*I_al+E_aln.*I_aln+E_si.*I_si;
rho_t = rho_aln.*t_aln + rho_si.*t_si + rho_al.*t_al;
N = (4.*E_tot*((16/3).*((m./2 + m_s).*L_b).^2 + 4.*(m.*L_m - 2.*m_s.*L_s).*(m./2 + m_s).*L_b + (m.*L_m - 2.*m_s.*L_s).^2));
D1 = (352/105).*(m./2 + m_s).^2.*(rho_t).*(w_b).*L_b.^5;
D2 = (104/45).*(m./2 + m_s).*((80/13).*m_s.^2+((80/13).*m-2.*rho_t.*w_b.*L_s).*m_s+m.*((20/13).*m+L_m.*w_b.*rho_t)).*L_b.^4;
D3 = (2/5).*(m.*L_m-2.*m_s.*L_s).*((160/3).*m_s.^2 + ((160/3).*m-2.*rho_t.*w_b.*L_s).*m_s+((40/3).*m+rho_t.*L_m.*w_b).*m).*L_b.^3;
D4 = 4.*(m./2 + m_s).*((40/3).*L_s.^2.*m_s.^2+m.*(L_m.^2-(34/3).*L_m.*L_s+L_s.^2).*m_s+(10/3).*L_m.^2.*m.^2).*L_b.^2;
D5 = 4.*(m.*L_m-2.*m_s.*L_s).*(4.*L_s.^2.*m_s.^2+m.*(L_m-L_s).^2.*m_s+L_m.^2.*m.^2).*L_b;
D6 = (L_m.^2.*m+2.*L_s.^2.*m_s).*(m.*L_m-2.*m_s.*L_s).^2;
D = L_b.*(D1+D2+D3+D4+D5+D6);
w2=N./D;
Freq =sqrt(w2)./(2.*pi)
if(Freq >= 810 || Freq <= 750)
obj = 10000000;
end
else
 obj = 100000;
end
end

[Sam Chak]

@Pavitra Jain, This is just a test to trigger the "If" condition. Note that there is another objective value "obj = Fom".

xTest = [500 5 100 1e8 10 100 100 5];
J     = objective(xTest)

w_b = 5.0000e-06

w_s = 50.0000

Fom = -0.0012

obj = -0.0012

Freq = 0.3471

J = 10000000

function obj = objective(x)
    % Defining Parameters
    L_b     = x(1)*10^-6;
    w_b     = x(2)*10^-6
    L_m     = x(3)*10^-6;
    w_m     = x(4)*10^-6;
    t_hSi   = x(5)*10^-6;
    t_aln   = x(6)*10^-9;
    L_s     = x(7)*10^-6;
    gap_s   = x(8)*10^-6;
    
    w_s     = 0.5*(w_m - w_b - 2*gap_s)
    t_si    = 1*10^-6;
    t_sio2  = 3*10^-6;
    t_al    = 0.1*10^-6;
    
    if (w_b <= 0.3.*L_b && w_s >= 5)
        %% Calculating Required parameters
        % Densities
        rho_al   = 2700;
        rho_aln  = 3300;
        rho_si   = 2329;
        rho_sio2 = 2200;
        
        % Elastic Modulus of Materials
        E_al     = 70*10^9;
        E_aln    = 330*10^9;
        E_si     = 188*10^9;
        
        % neutral axis
        n        = (E_al.*t_al.^2 + 2.*E_al.*t_al.*t_aln + 2.*E_al.*t_al.*t_si + E_aln.*t_aln.^2 + 2.*E_aln.*t_aln.*t_si +E_si.*t_si.^2)/(2.*(E_al.*t_al + E_aln.*t_aln + E_si.*t_si));
        I_al     = (1/12).*w_b.*t_al.^3 + w_b.*t_al.*(t_aln + t_si - n + 0.5.*t_al)^2;
        I_aln    = (1/12).*w_b.*t_aln.^3 + w_b.*t_aln.*(t_si - n + 0.5.*t_aln)^2;
        I_si     = (1/12).*w_b.*t_si.^3 + w_b.*t_si.*(0.5.*t_si - n)^2;
        m        = L_m.*w_m.*(rho_aln.*t_aln + rho_si.*t_si + rho_sio2.*t_sio2 + rho_si.*t_hSi);
        Force    = m.*9.8;
        d_31     = -1.9159*10^-12;
        m_s      = L_s.*w_s.*(rho_aln.*t_aln + rho_si.*t_si + rho_sio2.*t_sio2 + rho_si.*t_hSi);
        
        %% Charge Calculation
        Charge   = (d_31.*9.8.*w_b.*L_b.*t_aln./I_aln).*((0.5.*m+m_s)*L_b - m_s.*L_s + 0.5.*m.*L_m);
        
        %% Area Calculation
        Area     = (L_b+L_m).*w_m;
        
        %% Figure of Merit
        Normalised_Charge = Charge.*10^13;
        Normalised_Area   = Area.*10^9;
        Fom               = Normalised_Charge./Normalised_Area
        obj               = Fom
        
        %% Frequency
        E_tot   = E_al.*I_al+E_aln.*I_aln+E_si.*I_si;
        rho_t   = rho_aln.*t_aln + rho_si.*t_si + rho_al.*t_al;
        N       = (4.*E_tot*((16/3).*((m./2 + m_s).*L_b).^2 + 4.*(m.*L_m - 2.*m_s.*L_s).*(m./2 + m_s).*L_b + (m.*L_m - 2.*m_s.*L_s).^2));
        D1      = (352/105).*(m./2 + m_s).^2.*(rho_t).*(w_b).*L_b.^5;
        D2      = (104/45).*(m./2 + m_s).*((80/13).*m_s.^2+((80/13).*m-2.*rho_t.*w_b.*L_s).*m_s+m.*((20/13).*m+L_m.*w_b.*rho_t)).*L_b.^4;
        D3      = (2/5).*(m.*L_m-2.*m_s.*L_s).*((160/3).*m_s.^2 + ((160/3).*m-2.*rho_t.*w_b.*L_s).*m_s+((40/3).*m+rho_t.*L_m.*w_b).*m).*L_b.^3;
        D4      = 4.*(m./2 + m_s).*((40/3).*L_s.^2.*m_s.^2+m.*(L_m.^2-(34/3).*L_m.*L_s+L_s.^2).*m_s+(10/3).*L_m.^2.*m.^2).*L_b.^2;
        D5      = 4.*(m.*L_m-2.*m_s.*L_s).*(4.*L_s.^2.*m_s.^2+m.*(L_m-L_s).^2.*m_s+L_m.^2.*m.^2).*L_b;
        D6      = (L_m.^2.*m+2.*L_s.^2.*m_s).*(m.*L_m-2.*m_s.*L_s).^2;
        D       = L_b.*(D1+D2+D3+D4+D5+D6);
        w2      = N./D;
        Freq    = sqrt(w2)./(2*pi)        
        if (Freq >= 810 || Freq <= 750)
            obj = 10000000;
            % asd = 1;
        end        
    else
        obj = 100000;
        % asd = 1;
    end
end

Answer 1

[Torsten]

As you can see from the result, ga never seems to execute the "if" condition in your objective.

I don't know why you test for the constraint condition in the objective since you tried to implemented it in A and b in the call to "ga" ( but not correctly ).

[Pavitra Jain]

My constraint was meant to be x(1).*(-0.3) + x(2) <=0 and 2*x(8) + x(2) - x(4) <= -5. Thats what I have tried to specify in my If condition as well.

My problem here is when I run the code and when the code is executing,

ga(@objective,nvars,A,b,Aeq,beq,lb,ub,[],opts)

at this line for some reason when I put a breakpoint at the first line in the objective function, it never pauses indicating that it is not iterating the objective function for some reason.

If condition are meant as a fail safe in my earlier attempts I have seen that even though I have provided linear inequality but even then it still uses some of those sets of values which falls outside the constraints. These sets of value can result in NaN values as the equations are only valid for the linear inequality condition.

If condition for Freq also works as a constraint. I am trying to maximise Fom = Charge/Area for Freq value being in a certain range.

My main issue is the fact that it appears to me the objective function is not being executed at all. I am not concerned with the If condition. But Objective function not being executed is my real concern.

[Pavitra Jain]

One more thing @Torsten can you specify what mistake I have made in implementing in if condition and A and b for inequality.

[Torsten]

My constraint was meant to be x(1).*(-0.3) + x(2) <=0 and 2*x(8) + x(2) - x(4) <= -5. Thats what I have tried to specify in my If condition as well.

In your code you set the constraints

x(2)*10^-6 -0.3*x(1)*10^-6 <= 0
0.5*(x(4)*10^-6 - x(2)*10^-6 - 2*x(8)*10^-6) >=5

So the first inequality in the A and b representation is correct, the second not since b(2) should be -10*1e6 instead of -5.

Put all constraints in the A and b and remove the if-else constructs in the objective function.

[Pavitra Jain]

Ohh yes you are right I totally missed that.

[Pavitra Jain]

Thanks @Torsten but I think it still didnt resolve the main problem I am facing that is it never paused at the breakpoint I specified at the first line of the objective function.

[Pavitra Jain]

And yes it turns on the parallel pool but does not use it and leaves it idle.

[Torsten]

From my intuition, specifying an OutputFcn together with Parallel Computing makes no sense. Same for setting breakpoints in Parallel Computing. Setting breakpoints is debugging - the opposite of Parallel Computing. What do you think ?

I suggest you make your code work properly without Parallel Computing first and finally add this option as a cream topping.

[Pavitra Jain]

I have actually done this before to specify OutputFcn along with Parallel Computing. I always have the output function store the result at the end of each generation in a text file. I have never encountered such issue before. Same goes for debugging. Even if I consider your proposal for Breakpoint, parallel pool shouldn't go idle. Even if I dont use parallel pool, each run for objective function takes less than 30 ms. So, even without parallel pool, it would take 30s for 1000 iterations. So, it should store results of first generation in 30 secs time in the Output_Data.txt file.

[Torsten]

I always have the output function store the result at the end of each generation in a text file. I have never encountered such issue before.

Then you should think about what you've changed in your code that made it from working before to not working now. I cannot imagine that OutputFcn and Debugging can be combined with Parallel Computing - at least it shouldn't be done in my opinion.