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!
Best Answer