In this tutorial, you will learn how to use Matlab fminconfunction 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. Hint: double click the code snippets to select all. Note: the line numbers in each step is refer to as the code snippets in each step instead of the top3d program.
[expand title="Step.0: Remove Codes from Top3d" expanded="true" tag="h5"]

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.

[/expand]

 

[expand title="Step.1: Initialize Fmincon" expanded="true" tag="h5"]

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 both myObjFcn and myHessianFcn.
  • Since our constraint (volume constraint) is a function of xPhys instead of design variable x, we will define it in the function myConstrFcn
  • There are lots of OPTIONS defined for fmincon
    • TolX: user-defined terminarion criterion.
    • MaxIter: maximum number of iterations.
    • Algorithm: Matlab fmincon has many alogrithms, such as 'sqp', 'active-set', 'trust-region-reflective' and 'interior-point'. The reason why we choose 'interior-point' instead of others is because 'interior-point' accepted user-supplied Hessian of the Lagrange function while 'sqp' and 'active-set' do not allow user to provide hessian. The hessian in 'sqp' or 'active-set' are approximated by the means of BFGS (Quasi-Newtown 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 the OutputFcn.
    • OutputFcn: we defined output function in myOutputFcn, which will display iteration information, structural topology corresponding.
    • PlotFcns: Matlab provides some built-in plot function, @optimplotfval plots the function value.
    • Click here for more information out fmincon options.

[matlab title="Initialize Fmincon Matlab Code"]
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','interior-point',...
'GradObj','on', 'GradConstr','on', 'Hessian','user-supplied', 'HessFcn',@myHessianFcn,...
'Display','none', 'OutputFcn',@(x,optimValues,state) myOutputFcn(x,optimValues,state,displayflag), 'PlotFcns',@optimplotfval);
[/matlab]

[/expand]

[expand title="Step.2: Define Objective function" expanded="true" tag="h5"]

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

eq1

and the gradient is given by

eq

Click here to learn more about problem statement, objective function and sensitivity analysis.

[matlab light="false" title="Objective function: myObjFcn"]
function [f, gradf] = myObjFcn(x)
xPhys(:) = (Hx(:))./Hs;
% FE-ANALYSIS
sK = reshape(KE(:)
(Emin+xPhys(:)'.^penal*(E0-Emin)),2424nele,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
(E0-Emin)).ce)));
dc = -penal
(E0-Emin)*xPhys.^(penal-1).ce;
% FILTERING AND MODIFICATION OF SENSITIVITIES
dc(:) = H
(dc(:)./Hs);
% RETURN
f = c;
gradf = dc(:);
end % myfun
[/matlab]

See Writing Objective Functions for details.

[/expand]

[expand title="Step.3: Define Hessian" expanded="true" tag="h5"]

The Hessian of the Lagrangian is given by the equation:

Hessian of the Lagrangian 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

eq

Details about the derivation of Hessian corresponding to the objective function is discussed here.

myHessianFcn computes the Hessian at a point x with Lagrange multiplier structure lambda:

[matlab light="false" title="Hessian Function: myHessianFcn"]
function h = myHessianFcn(x, lambda)
xPhys = reshape(x,nely,nelx,nelz);
% Compute Hessian of Obj.
Hessf = 2*(penal*(E0-Emin)xPhys.^(penal-1)).^2 ./ (E0 + (E0-Emin)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
[/matlab]

For details, see Hessian. For an example, see fmincon Interior-Point Algorithm with Analytic Hessian.

[/expand]

[expand title="Step.4: Define Constraint" expanded="true" tag="h5"]

The function that computes the nonlinear inequality constraints c(x)≤ 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

eq

Click here to learn more about problem statement and constraints.

[matlab light="false" title="Constraint Function: myConstrFcn"]
function [cneq, ceq, gradc, gradceq] = myConstrFcn(x)
xPhys(:) = (Hx(:))./Hs;
% Non-linear Constraints
cneq = sum(xPhys(:)) - volfrac
nele;
gradc = ones(nele,1);
% Linear Constraints
ceq = [];
gradceq = [];
end % mycon
[/matlab]

For more information, see Nonlinear Constraints.
[/expand]

[expand title="Step.5: Define Output Function" expanded="true" tag="h5"]

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.

[matlab title="Output function: myOutputFcn"]
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
[/matlab]

See Output Function.

[/expand]

[expand title="Step.6: Call Fmincon" expanded="true" tag="h5"]

[matlab light="false" title="Call Fmincon Matlab Code"]fmincon(@myObjFcn, x, A, B, Aeq, Beq, LB, UB, @myConstrFcn, OPTIONS);[/matlab]

Click here to learn more about Matlab fmincon or excute the following command line in the Matlab:

[matlab light="true"]doc fmincon[/matlab]

[/expand]

[expand title="Step.7: Run the program" expanded="true" tag="h5"]
Now you can run the program, e.g., with the following Matlab command line:

[matlab light="true"]top3d(30,10,2,0.3,3,1.2)[/matlab]

[/expand]

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 01760-2098, UNITED STATES, 508-647-7000, Fax: 508-647-7001, info@mathworks.com, www.mathworks.com.