Speeding up scattered interpolation and my code in general

Hello Everyone,

I am using a genetic alorithm to run a series of beam dimension optimisation for various lengths. My issue is that my code seems to me to be very inefficient as it takes roughly 30 mins per beam length. I have to do a total of 150 different beam lengths so would take about 3 days to complete.

I think the process that takes up the most time is the scattered interpolant used in my nlcon function. It is called up at each iterations so becomes very time consuming. However, I am completely new to Matlab so struggling to think of a way to speed the scattered interpolation. Does anyone have an example of how to speed up a scattered interpolant within a loop?

Any other ways to speed up my code would be appreciated as well. I have shown my code below doing 1 iteration of beam length (L) of 2300.

Thank you for your time and hope my code isn't too messy to read!

main script:

clc,clear all, close all
nvars = 4;                                 %Initial guess
A =[];
b = [];
Aeq = [];
beq = [];
lb = [50; 5; 100; 4];
ub = [450; 70; 1100; 36];

L0 = 2300:100:2300;
x = zeros(4,1);
for i=1:numel(L0)
L =L0(i);
%options for the analysis
hybridopts = optimoptions('fmincon','Display','iter',...
    'MaxFunctionEvaluations', 1000000, 'MaxIterations',10000,...
    'FunctionTolerance',1e-6,'OptimalityTolerance',1e-6);
options =optimoptions('ga','Display','iter','ConstraintTolerance',1e-6,...
    'FunctionTolerance',1e-6, 'PopulationSize',150,'PlotFcn', @gaplotbestf,...
    'HybridFcn',{'fmincon',hybridopts});  
%optimisation
[x(:,i)] = ga(@(x)objective(x), nvars, A,b,Aeq,beq,lb,ub,@(x)nlcon(x,L), options);
end

Objective function:

function f = objective(x) 
f= ((2*x(1)*x(2)) + ((x(3)-(2*x(2)))*x(4))); %Minimise cross sectional area

nlcon function:

function [c,ceq] = nlcon(x,L)

% Initial information
w = 51.1;                               % N/mm load
                                        % mm Span of beam
E = 210000;                             % N/mm2 Youngs Modulus
v = 0.3;                                % Poisson's Ratio  
ymo = 1;                                % Material Partial Factor

% Yield strength
if max([x(2),x(4)]) <= 16               % N/mm2 yield strength for S355 JR/J0 steel
    fy =355;
elseif max([x(2),x(4)]) <= 40
    fy =345;
elseif max([x(2),x(4)]) <= 63
    fy = 335;
elseif max([x(2),x(4)]) <= 80
    fy = 325;
elseif max([x(2),x(4)]) <= 100
    fy = 315;
else
    fy = 295;
end

% Plastic section modulus
A1 = x(1)*x(2);                         % mm2 area of single flange 
A2 = ((x(3)-(2*x(2)))/2)*x(4);          % mm2 area of half a flange
y1 = ((x(3)-(2*x(2)))/2)+(x(2)/2);      % mm from PNA to centre of the flange
y2 = ((x(3)-(2*x(2)))/4);               % mm from PNA to centre of the web
yp = ((A1*y1)+(A2*y2))/(A1+A2);         % mm centre of compression area
zp = ((A1+A2)*yp)+((A1+A2)*yp);         % mm3 Plastic section modulus and as compresson and tension flange are the same it is eseentially doubling it.

% Second Moment of area
I1 = ((x(1)*x(2)^3)/12)+((x(1)*x(2))*((((x(3)-(2*x(2)))/2)+(x(2)/2))^2)); % mm4 I+Ay2 Parallel axis thereom
I2 = (x(4)*(x(3)-(2*x(2)))^3)/12;                                         % I = bh3/12             
I3 = ((x(1)*x(2)^3)/12)+((x(1)*x(2))*((((x(3)-(2*x(2)))/2)+(x(2)/2))^2)); % same as I1 but for bottom flange
Ix = I1+I2+I3;                          % mm4 second moment of area of the I beam 

% For section classification
eps= (235/fy)^(1/2);                    % required for section classification

% Shear area of I-beam when load is parallel to web
Av = x(3)*x(4);                         % mm2 the web area when load is parallel to web

% Section Classification Flange width to thickness ratio in compression
% Class 2 taken to be the upper limit as same procedure as Class 1 
c1 = (((x(1)-x(4))/2)/(x(2)))-(10*eps);

% Section Classification: Web width to thickness ratio in bending
c2 = (x(3)-(2*(x(2))))/(x(4))-(83*eps);

% Bending Moment Constraint
Emx = (w*1.5)*L^2/2;                     % Nmm Bending Moment of a UDL Cantilever Beam wl2/2
Rmx = (zp*fy)/ymo;                       % Nmm Bending Resistance in accordance to Eurocode 3 for Class 1 or 2 beams 
c3 = Emx-Rmx;                            % Inequality constraint for Bending Moment

% Shear Force Constraint
Evx = (w*1.5)*L;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           % N Shear force of a UDL Cantilever beam  
Rvx = (Av*(fy*(3^(1/2))))/ymo;           % N Shear Resistance in accordance to Eurocode 3 for class 1 or 2
c4  = Evx-Rvx;

% Deflection Constraint
Edf = ((w*L^4)/(8*E*Ix));
Rdf = L/250;
c5 = Edf-Rdf;

%% Lateral-Torsional Buckling
% Section Properties
G = E/((2*(1+v)));                       % N/mm2 Shear Modulus in accordance with Eurocode 3 cl 3.2.6
It = (1/3)*(((x(3)-(2*x(2)/2))*(x(4)^3))+(2*x(1)*(x(2)^3))); % Torsional Constant in accordance to midas civil section manual

%Minor axis Second moment of area
Iz1 = (x(2)*x(1)^3)/12;
Iz2 = ((x(3)-(2*x(2)))*x(4)^3)/12; 
Iz3 = (x(2)*x(1)^3)/12;
Iz = Iz1 + Iz2 + Iz3;

%Warping Constant
Iw = (Iz*((x(3)-x(2))^2))/4;             %mm6 from SCI P385 cl B.3.1

%Warping rigidity of beam
kwt0 = (1/L)*(((E*Iw)/(G*It))^(1/2));

% n parameter showing destablising or stablising load applied.
za = x(3)/2;                             % assumed load applied on top fibre of beam
hs = x(3)-x(2);                          % distance between shear centre of flange
n0 = za/(hs/2);

% Calculating C using P385 SCI

n=[-2;-2;-2;-2;-2;-2;-2;-2;-2;-2;-1.5;-1.5;-1.5;-1.5;-1.5;-1.5;-1.5;-1.5;-1.5;-1.5;...
    -1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-0.5;-0.5;-0.5;-0.5;-0.5;-0.5;-0.5;-0.5;-0.5;-0.5;...
    0;0;0;0;0;0;0;0;0;0;0.25;0.25;0.25;0.25;0.25;0.25;0.25;0.25;0.25;0.25;...
    0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.75;0.75;0.75;0.75;0.75;0.75;0.75;0.75;0.75;0.75;...
    1;1;1;1;1;1;1;1;1;1;1.25;1.25;1.25;1.25;1.25;1.25;1.25;1.25;1.25;1.25;...
    1.5;1.5;1.5;1.5;1.5;1.5;1.5;1.5;1.5;1.5;2;2;2;2;2;2;2;2;2;2;...
    2.5;2.5;2.5;2.5;2.5;2.5;2.5;2.5;2.5;2.5;3;3;3;3;3;3;3;3;3;3];
kwt=[0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;...
    0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;...
    0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;...
    0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;...
    0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;...
    0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;...
    0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1;0;0.05;0.1;0.15;0.2;0.3;0.4;0.6;0.8;1];
c=[2.04;2.42;2.87;3.37;3.93;5.13;6.4;9.07;11.8;14.6;2.04;2.34;2.71;3.14;3.62;...
    4.67;5.79;8.16;10.6;13.1;2.04;2.25;2.53;2.87;3.25;4.11;5.04;7.04;9.11;11.2;...
    2.04;2.16;2.34;2.56;2.82;3.39;4.02;5.38;6.82;8.3;2.04;2.06;2.13;2.22;2.32;...
    2.5;2.66;2.88;3;3.08;2.04;2.02;2.03;2.05;2.06;2.06;2.02;1.85;1.65;1.47;...
    2.04;1.97;1.92;1.87;1.82;1.69;1.54;1.26;1.03;0.87;2.04;1.92;1.82;1.71;1.6;...
    1.39;1.21;0.92;0.73;0.61;2.04;1.87;1.71;1.56;1.41;1.17;0.98;0.72;0.56;0.46;...
    2.04;1.82;1.61;1.42;1.25;0.99;0.81;0.59;0.46;0.37;2.04;1.77;1.52;1.3;1.12;...
    0.86;0.7;0.5;0.38;0.31;2.04;1.67;1.35;1.09;0.91;0.68;0.54;0.38;0.29;0.23;...
    2.04;1.58;1.2;0.93;0.76;0.55;0.43;0.3;0.23;0.19;2.04;1.49;1.07;0.81;0.65;...
    0.47;0.36;0.25;0.19;0.16];
F = scatteredInterpolant(n, kwt, c);
C1 = max([F(n0,kwt0),0.001]);

%Critical Buckling moment for Fork supported simply supported beam
Mcr0 = ((pi)/L)*((E*Iz*G*It)^(1/2));

%Critical buckling Load for cantilever UDL beam.
Mcr = C1*Mcr0;

% Imperfection factor alphaLT
if x(3)/x(1)<=2
    alphaLT = 0.49;
else
    alphaLT = 0.79;
end


%Non-dimensional Slenderness lambdaLT calculation
lambdaLT = ((zp*fy)/Mcr)^(1/2);
lambdaLT0 = 0.2;                   % From NA to BS EN 1993-1-1 cl. NA.2.17b for welded sections
BetaLT = 1;                        % From NA to BS EN 1993-1-1 cl. NA.2.17b for welded sections

if (lambdaLT<=lambdaLT0) || Emx/Mcr <= (lambdaLT0^2)
    c6 = Emx-Rmx;
else
%phiLT
phiLT = 0.5*(1+(alphaLT*(lambdaLT-lambdaLT0))+(BetaLT*lambdaLT^2));

%Reduction Factor 

XLT0 = 1/(phiLT+(((phiLT^2)-(BetaLT*(phiLT^2)))^(1/2)));
XLT = min([XLT0,1]);
%Buckling Resistance
Mbrd = XLT*zp*fy/ymo;

c6 = Emx-Mbrd;
end

c7 = (x(3)/x(1))-3.1;

if Evx/Rvx > 0.5
    p = (((2*Evx)/Rvx)-1)^2;
    fyr = (1-p)*fy;
    Rmxr = (zp*fyr)/ymo;
    c8 = Emx-Rmxr;
else c8 = Emx-Rmx;
end
%% contraint inequality
c = [c1, c2, c3, c4, c5, c6, c7, c8];

%% constraint equality
ceq=[];

[darova]

You definitely should create your interpolant outside function

% data
F = scatterredInterpolant(data);
[x(:,i)] = ga(@(x)objective(x), nvars, A,b,Aeq,beq,lb,ub,@(x)nlcon(x,L,F), options);