In this tutorial, you will learn how to use Matlab1
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.
Note: the line numbers in each step is refer to as the code snippets in each step instead of the
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.
In this step, we are going to initialize parameters used by Matlab
- In line 2, a global variable
cewill be used in both
- Since our constraint (volume constraint) is a function of
xPhysinstead of design variable
x, we will define it in the function
- There are lots of
TolX: user-defined terminarion criterion.
MaxIter: maximum number of iterations.
fminconhas 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).
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: we defined output function in
myOutputFcn, which will display iteration information, structural topology corresponding.
PlotFcns: Matlab provides some built-in plot function,
@optimplotfvalplots the function value.
- Click here for more information out
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
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 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 is a vector that contains the nonlinear inequalities evaluated at
ceq is a vector that contains the nonlinear equalities evaluated at
gradc, the gradient of
gradceq, the gradient of
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
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.
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:
Step.7: Run the program
Now you can run the program, e.g., with the following Matlab command line:
If you have any questions, difficulties with this tutorial, please don’t hesitate to contact us.