MATLAB Answers

using fmincon with extra parameters

2 views (last 30 days)
Sadu on 31 Mar 2020
Edited: Sadu on 1 Apr 2020
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);
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;
alphaLT = 0.79;
%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
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
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!


Sign in to comment.

Accepted Answer

Cris LaPierre
Cris LaPierre on 31 Mar 2020
See this post from the MathWorks Support Team.
You'd probably need a loop to iterate through all the lengths. Update your call to fmincon like so
L = 4300;
x = fmincon(objective, x0, A,b,Aeq,beq,lb,ub,@(x)nlcon(x,L),options);
function [c,ceq] = nlcon(x,L)


Show 2 older comments
Sadu on 1 Apr 2020
I have changed the objective as I wanted the objective function to be seperate to the main code.
But I am running into the same problem. On the 11th loop i get the same error. I am going to start having a look at it to see if i can fix it but I do not understand why it can do 10 loops correctly but at the 11th it fails to make the constraint function work.
The first 10 values are as expected so it is looking good so far!
Sadu on 1 Apr 2020
I have not been successful in finding a solution to this problem. I have identified the error occurs for any length (L) above 1020 but other than that I do not know what I could do to fix this. I do not believe values above 1020 would give a complex number so unsure how the error occurs.
Edit: I have run the non linear constraint code seperately using the initial values and L=1100 with no issue getting the constraint values so no idea why I get the error mentioned above
Sadu on 1 Apr 2020
I have located the issue, which was kwt0 is large when the fmincon iterates at the smaller dimensions so the C1 value returned from the scattered interpolant is negative, which makes the Mbrd a complex number.
Is there anyway to make fmincon skip the optimisation iteration with ktw0 values >1?
Note: i do not mean skip the for loop iteration for the length but the actual fmincon iteration so it continues finding the optimal value while skipping all iterations with kwt0 >1.
Edit: I have sorted this out now and my code works, though I think I have to move to genetic algorithm as the solutions seem to be the local minimum around my initial guess, which gives non-optimal values when length reaches past 10000 so thought a genetic algorithm may work better in this case.
Thank you for taking time to help me and I hope you have a nice day!

Sign in to comment.

More Answers (0)

Sign in to answer this question.