using fmincon with extra parameters

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!

Answer 1

[Cris LaPierre]

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; 

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

...

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

[Sadu]

Cris, thank you for your reply, I have successfully passed the additional parameter (length) but now I am facing a 'In an assignment A(:) = B, the number of elements in A and B must be the

same.' error for the x(i) = fmincon(..) line.

I do understand this arised due to there being both scalar and arrays but I do not know how to overcome this issue. Could I please get any advice or get pointed towards a similar example please so i can amend my code.

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

L0 = 100:100:10000;
x = zeros(4,100);
for i=1:numel(L0)
L =L0(i);
%options for the analysis
options = optimoptions('fmincon', 'Algorithm', 'sqp', 'Display', 'iter-detailed',...
    'MaxFunctionEvaluations', 1000000, 'MaxIterations', 20000,...
    'FunctionTolerance', 1e-1);

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

[Cris LaPierre]

The issue is with your assignment. You are trying to assing 4 values to a single element of x. Modify your assignment to include row and columns. Based on how you preallocate x, this should work.

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

Have you changed your objective function? If I use the objective function you originally shared, I don't need the @(x) and the (x,L) around objective.

Also, if you are going to run that many loops, you are going to want to turn off the printing to screen (add semicolon to the end of lines that are printing values to screen). Writing to screen slows down your code.

I hit a new error in the 11th loop.

Error using sqpInterface
Nonlinear constraint function is undefined at initial point. Fmincon cannot continue.

Error in fmincon (line 843)
    [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = sqpInterface(funfcn,X,full(A),full(B),full(Aeq),full(Beq), ...

[Sadu]

Cris,

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]

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]

Cris,

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!