3D Topology Optimization using Method of Moving Asymptotes - Top3d/MMA

In this tutorial, you will learn how to implement Method of Moving Asymptotes(MMA) to our efficient 3d topology optimization program (top3d) in Matlab¹.

Updates:

• The parameter c_MMA is increased by a power of 10.  The gradient of the mass constraint dfdx is divided by (volfrac * nele).
• The MMA solver used in this tutorial is Version September 2007. (Aug 2017)

Step.1: Obtain MMA Matlab files

Before you can use Method of Moving Asymptotes (MMA) as an optimizer in our 3d topology optimization program, you need to obtain the Matlab implementation of MMA from Prof. Krister Svanberg (krille@math.kth.se) from KTH in Stockholm Sweden.

Step.2: Initialize MMA

Suppose now you have the MMA algorithm, which contains mmasub.m and subsolv.m. We can start modifying our program.

There is total of 29 input and output variables to mmasub function. Some of them can be predefined before entering the iteration loop. The following code can be interred after line 65 in the top3d program. For details about how to chose constants like a0, a, c_MMA, d, please check the Matlab note along with the MMA Matlab algorithm.

% INITIALIZE MMA OPTIMIZER
m     = 1;                % The number of general constraints.
n     = nele;             % The number of design variables x_j.
xmin  = zeros(n,1);       % Column vector with the lower bounds for the variables x_j.
xmax  = ones(n,1);        % Column vector with the upper bounds for the variables x_j.
xold1 = x(:);             % xval, one iteration ago (provided that iter>1).
xold2 = x(:);             % xval, two iterations ago (provided that iter>2).
low   = ones(n,1);        % Column vector with the lower asymptotes from the previous iteration (provided that iter>1).
upp   = ones(n,1);        % Column vector with the upper asymptotes from the previous iteration (provided that iter>1).
a0    = 1;                % The constants a_0 in the term a_0*z.
a     = zeros(m,1);       % Column vector with the constants a_i in the terms a_i*z.
c_MMA = 10000*ones(m,1);  % Column vector with the constants c_i in the terms c_i*y_i.
d     = zeros(m,1);       % Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.

Step.3: Implement MMA

The implementation of MMA algorithm is very simple, one just need to replace the OC algorithm in lines 81-88 with MMA algorithm. However, there are some auxiliary variables we need define before we call the mmasub function. Those are defined in the following code:

% METHOD OF MOVING ASYMPTOTES
xval  = x(:);
f0val = c;
df0dx = dc(:);
fval  = sum(xPhys(:))/(volfrac*nele) - 1;
dfdx  = dv(:)' / (volfrac*nele);
[xmma, ~, ~, ~, ~, ~, ~, ~, ~, low,upp] = ...
mmasub(m, n, loop, xval, xmin, xmax, xold1, xold2, ...
f0val,df0dx,fval,dfdx,low,upp,a0,a,c_MMA,d);
% Update MMA Variables
xnew     = reshape(xmma,nely,nelx,nelz);
xPhys(:) = (H*xnew(:))./Hs;
xold2    = xold1(:);
xold1    = x(:);

Please note that, in the code shown above, the inequality constraint is normalized. Also, to achieve better results, please consider smaller value than default tolx.

FAQ

Why my volume fraction is not satisfied?

By following this tutorial, the MMA is implemented with inequality constraint (as defined in mmasub.m), which is in the form of $Ax-b\le 0$. In order to impose equality constraint (e.g., $Ax-b=0$), you will need two constraints: i.e., $Ax-b\le 0$, $-Ax+b\ge 0$.

If you have any questions, difficulties with this tutorial, please don’t hesitate to contact us.