3D Topology Optimization using MATLAB fmincon – Top3d/fmincon

In this tutorial, you will learn how to use Matlab¹ 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 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.

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);

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

$c(\tilde{\mathbf{x}}) = \mathbf{F}^{\text{T}}\mathbf{U}(\tilde{\mathbf{x}}),$

and the gradient is given by

$\frac {\partial c(\tilde{\mathbf{x}})} {\partial \tilde x_{i}} = - \mathbf u_{i}(\tilde{\mathbf{x}})^{\text{T}} \left [p \tilde x_{i}^{p-1} (\text{E}_{0} - \text{E}_{\min}) \mathbf{k}_{i}^{0}\right ] \mathbf u_{i}(\tilde{\mathbf{x}}).$

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

function [f, gradf] = myObjFcn(x)
    xPhys(:) = (H*x(:))./Hs;
    % FE-ANALYSIS
    sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(E0-Emin)),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*(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

See Writing Objective Functions for details.

Step.3: Define Hessian

The Hessian of the Lagrangian is given by the equation:

$\nabla_{xx}^2L(x, \lambda) = \nabla^2 f(x) + \sum \lambda_i \nabla^2 c_i(x) + \sum \lambda_i \nabla^2 ceq_i(x).$

In our problem the Hessian of the objective function can be expressed as:

since the constraint is linear constraint, the Hessian of the constraint

$\nabla^2 c_i (\tilde{\mathbf{x}}) = 0$

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*(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

For details, see Hessian. For an example, see fmincon Interior-Point 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

$v(\tilde{\mathbf{x}}) = \tilde{\mathbf{x}}^{\mathrm{T}} \mathbf{v} - \bar{v} \leqslant 0$

Click here to learn more about problem statement and constraints.

function [cneq, ceq, gradc, gradceq] = myConstrFcn(x)
    xPhys(:) = (H*x(:))./Hs;
    % Non-linear 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 01760-2098, UNITED STATES, 508-647-7000, Fax: 508-647-7001, info@mathworks.com, www.mathworks.com. ↩