3D Topology Optimization using MATLAB fmincon – Top3d/fmincon
In this tutorial, you will learn how to use Matlab^{1}
fmincon
function as an optimizer in our 3d topology optimization program. It constrains six(6) main steps, i.e., Initialize Fmincon, Define Objective function, Hessian, Constraint, Output function and Call fmincon.
Note: the line numbers in each step is refer to as the code snippets in each step instead of the top3d
program.
Step.0: Remove Codes from Top3d
Before we get started, please download the Top3d
program if you haven’t done so.
Then, DELETE lines 64 96 (mainly the while
loop) in the program. The following steps start after line 63.
Step.1: Initialize fmincon
In this step, we are going to initialize parameters used by Matlab fmincon
function.
 In line 2, a global variable
ce
is defined.ce
will be used in bothmyObjFcn
andmyHessianFcn
.  Since our constraint (volume constraint) is a function of
xPhys
instead of design variablex
, we will define it in the functionmyConstrFcn
 There are lots of
OPTIONS
defined forfmincon
TolX
: userdefined terminarion criterion.MaxIter
: maximum number of iterations.Algorithm
: Matlabfmincon
has many alogrithms, such as ‘sqp’, ‘activeset’, ‘trustregionreflective’ and ‘interiorpoint’. The reason why we choose ‘interiorpoint’ instead of others is because ‘interiorpoint’ accepted usersupplied Hessian of the Lagrange function while ‘sqp’ and ‘activeset’ do not allow user to provide hessian. The hessian in ‘sqp’ or ‘activeset’ are approximated by the means of BFGS (QuasiNewtown method).GradObj,
GradConstr
,Hessian
,HessFcn
: With analytic expressions of the gradients of the objective function, constraints and Hessian of the Lagrange function, the optimization algorithm will execute faster than numerical approximation.Display
: display option is turned off. We will display iteration information in theOutputFcn
.OutputFcn
: we defined output function inmyOutputFcn
, which will display iteration information, structural topology corresponding.PlotFcns
: Matlab provides some builtin plot function,@optimplotfval
plots the function value. Click here for more information out
fmincon
options.
global ce % Shared between myfun and myHessianFcn
A = [];
B = [];
Aeq = [];
Beq = [];
LB = zeros(size(x));
UB = ones(size(x));
OPTIONS = optimset('TolX',tolx, 'MaxIter',maxloop, 'Algorithm','interiorpoint',...
'GradObj','on', 'GradConstr','on', 'Hessian','usersupplied', 'HessFcn',@myHessianFcn,...
'Display','none', 'OutputFcn',@(x,optimValues,state) myOutputFcn(x,optimValues,state,displayflag), 'PlotFcns',@optimplotfval);
Step.2: Define Objective function
The function to be minimized. myObjFcn
is a function that accepts a vector x
and returns a scalar f
and the gradient vector gradf(x)
, the objective function evaluated at x
.
In our problem (minimize compliance) the objective is given by
and the gradient is given by
Click here to learn more about problem statement, objective function and sensitivity analysis.
function [f, gradf] = myObjFcn(x)
xPhys(:) = (H*x(:))./Hs;
% FEANALYSIS
sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(E0Emin)),24*24*nele,1);
K = sparse(iK,jK,sK); K = (K+K')/2;
U(freedofs,:) = K(freedofs,freedofs)F(freedofs,:);
% OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
ce = reshape(sum((U(edofMat)*KE).*U(edofMat),2),[nely,nelx,nelz]);
c = sum(sum(sum((Emin+xPhys.^penal*(E0Emin)).*ce)));
dc = penal*(E0Emin)*xPhys.^(penal1).*ce;
% FILTERING AND MODIFICATION OF SENSITIVITIES
dc(:) = H*(dc(:)./Hs);
% RETURN
f = c;
gradf = dc(:);
end % myfun
See Writing Objective Functions for details.
Step.3: Define Hessian
The Hessian of the Lagrangian is given by the equation:
In our problem the Hessian of the objective function can be expressed as:
since the constraint is linear constraint, the Hessian of the constraint
Details about the derivation of Hessian corresponding to the objective function is discussed in the paper.
myHessianFcn
computes the Hessian at a point x with Lagrange multiplier structure lambda:
function h = myHessianFcn(x, lambda)
xPhys = reshape(x,nely,nelx,nelz);
% Compute Hessian of Obj.
Hessf = 2*(penal*(E0Emin)*xPhys.^(penal1)).^2 ./ (E0 + (E0Emin)*xPhys.^penal) .* ce;
Hessf(:) = H*(Hessf(:)./Hs);
% Compute Hessian of constraints
Hessc = 0; % Linear constraint
% Hessian of Lagrange
h = diag(Hessf(:)) + lambda.ineqnonlin*Hessc;
end % myHessianFcn
For details, see Hessian. For an example, see fmincon
InteriorPoint Algorithm with Analytic Hessian.
Step.4: Define Constraint
The function that computes the nonlinear inequality constraints $c(x) \leq 0$ and the nonlinear equality constraints $ceq(x) = 0$. myConstrFcn
accepts a vector x
and returns the four vectors c
, ceq
, gradc
, gradceq
. c
is a vector that contains the nonlinear inequalities evaluated at x
, ceq
is a vector that contains the nonlinear equalities evaluated at x
, gradc
, the gradient of c(x)
, and gradceq
, the gradient of ceq(x)
.
In our problem, we only have one equality constraint, i.e., volume fraction constraint
Click here to learn more about problem statement and constraints.
function [cneq, ceq, gradc, gradceq] = myConstrFcn(x)
xPhys(:) = (H*x(:))./Hs;
% Nonlinear Constraints
cneq = sum(xPhys(:))  volfrac*nele;
gradc = ones(nele,1);
% Linear Constraints
ceq = [];
gradceq = [];
end % mycon
For more information, see Nonlinear Constraints.
Step.5: Define Output Function
The myOutputFcn
is a function that an optimization function calls at each iteration. The iteration information is displayed and the structural topology can be also displayed during the iteration process depends on users setting.
function stop = myOutputFcn(x,optimValues,state,displayflag)
stop = false;
switch state
case 'iter'
% Make updates to plot or guis as needed
xPhys = reshape(x, nely, nelx, nelz);
%% PRINT RESULTS
fprintf(' It.:%5i Obj.:%11.4f Vol.:%7.3f ch.:%7.3f\n',optimValues.iteration,optimValues.fval, ...
mean(xPhys(:)),optimValues.stepsize);
%% PLOT DENSITIES
if displayflag, figure(10); clf; display_3D(xPhys); end
title([' It.:',sprintf('%5i',optimValues.iteration),...
' Obj. = ',sprintf('%11.4f',optimValues.fval),...
' ch.:',sprintf('%7.3f',optimValues.stepsize)]);
case 'init'
% Setup for plots or guis
if displayflag
figure(10)
end
case 'done'
% Cleanup of plots, guis, or final plot
figure(10); clf; display_3D(xPhys);
otherwise
end % switch
end % myOutputFcn
See Output Function.
Step.6: Call fmincon
fmincon(@myObjFcn, x, A, B, Aeq, Beq, LB, UB, @myConstrFcn, OPTIONS);
Click here to learn more about Matlab fmincon
or execute the following command line in the Matlab:
doc fmincon
Step.7: Run the program
Now you can run the program, e.g., with the following Matlab command line:
top3d(30,10,2,0.3,3,1.2)
If you have any questions, difficulties with this tutorial, please don’t hesitate to contact us.

MATLAB is a registered trademark of The MathWorks, Inc. For MATLAB product information, please contact The MathWorks, Inc., 3 Apple Hill Drive, Natick, MA 017602098, UNITED STATES, 5086477000, Fax: 5086477001, info@mathworks.com, www.mathworks.com. ↩