Indiana University-Purdue University Indianapolis

3D Topology Optimization using MATLAB fmincon – top3d/fmincon

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.

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

eq1

and the gradient is given by

eq

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:

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:

    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)≤ 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.

    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 excute 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.

  • Was this Helpful ?
  • Yes   No