2 views (last 30 days)

I am in the process of writing an optimisation code for the design of beams. I am using the fmincon function to achieve the minimal cross sectional area. I have completed the code for 1 iteration, which is shown below:

The objective function:

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

x0 = [210,17.4,536.7,10.8]; %Initial guess

disp(['Initial Objective: ' num2str(objective(x0))]) %Initial guess C.A

%fmincon inputs

A =[];

b = [];

Aeq = [];

beq = [];

lb = [50; 5; 100; 4];

ub = [450; 70; 1100; 36];

%non linear inequalities and equalities constraints

nonlincon = @nlcon;

%options for the analysis

options = optimoptions('fmincon', 'Algorithm', 'sqp', 'Display', 'iter-detailed',...

'MaxFunctionEvaluations', 1000000, 'MaxIterations', 20000,...

'FunctionTolerance', 1e-1);

%optimisation

x = fmincon(objective, x0, A,b,Aeq,beq,lb,ub,nonlincon,options);

%Display final results with constraints shown

disp(['Final Objective: ' num2str(objective(x))])

[c,ceq] = nlcon(x)

the non-linear constraints:

function [c,ceq] = nlcon(x)

% Initial information

w = 51.1; % N/mm load

L = 4300; % mm Span of beam

E = 210000; % N/mm2 Youngs Modulus

v = 0.3; % Poisson's Ratio

ymo = 1; % Material Partial Factor

...

% 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;

...

...

...

...

% 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);

%phiLT

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

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;

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

%% contraint inequality

c = [c1, c2, c3, c4, c5, c6, c7];

%% constraint equality

ceq=[];

My question is how would I go about doing multiple iteration of this with different Length of beam (L)?

I want to find the optimal design variables for each length of beam between 100:100:15000. I was thinking about using a for loop but not sure how and where to implement this within the code. Any help or advice would be greatly appreciated. Thank you for your time!

Cris LaPierre
on 31 Mar 2020

You'd probably need a loop to iterate through all the lengths. Update your call to fmincon like so

...

L = 4300;

%optimisation

x = fmincon(objective, x0, A,b,Aeq,beq,lb,ub,@(x)nlcon(x,L),options);

...

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

...

end

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.