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

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; % 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:

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

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

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.