Optimality Criteria Method

Question

This MATLAB code shows how one can build a bridge according to topological optimization theory. The code is also based on Euler-Bernoulli beam theory.

The very last part (%% Optimality Criteria Update) consists of optimality criteria method. However, I have some trouble understanding this sentence.

xnew=max(0.001,max(x-move,min(1,min(x+move,x.*sqrt(-dc./lmid)))));

At the end of this sentence, one is going to compare 1, x+move, and x.*sqrt(-dc./lmid), and find the minimum value among them. However, I don't understand what this x.*sqrt(-dc./lmid) means, and why it suddenly comes out.

Can anyone explain what this x.*sqrt(-dc./lmid) means, and why do I have to compare it? This image might be of help.

%% Ground-structure based Topology Optimization by D.-M. Kim, May 2016
function gstop(Nx,Ny,Vlim,penal,maxL)
% Generate Nodes & Elements
[node,elem]=GS(Nx,Ny,maxL);

Ne=size(elem,1);    % Total number of elements
% Initialize
x=zeros(Ne,1);
x(1:Ne)=0.1;     % initial density
loop=0;
change=1;
change2=1;
c=0;
cold=0;
diffc=1;
diffc2=1;
xold=x;

filedir = 'images/Plot_iteration_';

% Start Iteration
while ((change>0.001 && change2>0.001) || (diffc>1e-2 && diffc2>1e-2))
    loop=loop+1;
    xold2=xold;
    xold=x;
    cold2=cold;
    cold=c;
    % FE-analysis
    [U]=FE(Nx,Ny,node,elem,x,penal);
    % Objective Function and Sensitivity Analysis
    [KE]=stiffness(elem);       % element stiffness matrix
    c=0.;
    dc=zeros(Ne,1);
    for jj=1:Ne
        n1=elem(jj,1);
        n2=elem(jj,2);
        Ue=U([3*n1-2;3*n1-1;3*n1;3*n2-2;3*n2-1;3*n2],1);
        c=c+x(jj)^penal*Ue'*KE(:,:,jj)*Ue;
        dc(jj)=-penal*x(jj)^(penal-1)*Ue'*KE(:,:,jj)*Ue;
    end
    % Design Update by the Optimality Criteria Method
    [x]=OC(x,elem,Vlim,dc);
    % Print Results
    change=max(abs(x-xold));
    change2=max(abs(x-xold2));
    diffc=abs((c-cold)/c);
    diffc2=abs((c-cold2)/c);
    disp(['It.:' sprintf('%4i',loop) ' / Obj.:' sprintf('%10.4f',c)...
        ' / Vol.:' sprintf('%6.3f',x'*elem(:,3))...
        ' / ch.:' sprintf('%6.3f',change)])
    % Plot Densities
    clf
    hold on
    plot(node(:,1),node(:,2),'.','color','r','markersize',3)
    [xsort,order]=sort(x);
    for jj=1:Ne
        if xsort(jj)>0.01
            n1=elem(order(jj),1);
            n2=elem(order(jj),2);            
            % undeformed
            x1=node(n1,1); y1=node(n1,2);
            x2=node(n2,1); y2=node(n2,2);
              % deformed
              x1=node(n1,1)+U(3*n1-2); y1=node(n1,2)+U(3*n1-1);
              x2=node(n2,1)+U(3*n2-2); y2=node(n2,2)+U(3*n2-1);
            r=3*xsort(jj);
            rr=1-xsort(jj);
            plot([x1,x2],[y1,y2],'linewidth',r,'color',rr*[0 0 0])

            % plotting the mirror image
            hold on
            plot([-x1, -x2], [y1, y2],'linewidth',r,'color',rr*[1 1 1])
        end
    end
    axis equal; axis tight; axis off; pause(1e-6);

    print(strcat(filedir, int2str(loop)), '-dpng');
end

%% Generate Ground-structure (Nodes & Elements)
function [node,elem]=GS(Nx,Ny,maxL)
Nn=Nx*Ny;       % Total Number of nodes
node=zeros(Nn,2);   % Location of nodes (col1: x-dir / col2: y-dir)
Lx=80; Ly=40;
% Generate Nodes
no=0;
for jy=1:Ny
    for jx=1:Nx
        no=no+1;
        node(no,1)=(Lx/(Nx-1))*(jx-1);    % x-dir. position
        node(no,2)=(Ly/(Ny-1))*(jy-1);    % y-dir. position
    end
end

% Generate Elements
elem=zeros(10000,4);    % Element data (node 1 / node 2 / elem. length / rotation angle)
no=0;
for n1=1:Nn
    no_temp=0;
    the_temp=[];
    for n2=n1+1:Nn
        x1=node(n1,1); y1=node(n1,2);
        x2=node(n2,1); y2=node(n2,2);
        dx=x2-x1; dy=y2-y1;
        E_len=norm([dx,dy]);
        E_theta=acos(dx/E_len);
        max_length=maxL*(Lx/(Nx-1));
        % Save element data only if the length is shorter than maxL
        if (E_len<=max_length)&&(no_temp==0||isempty(find(abs(the_temp-E_theta)<1e-3,1)))
            no=no+1;
            elem(no,:)=[n1,n2,E_len,E_theta];
            no_temp=no_temp+1;
            the_temp(no_temp,1)=E_theta;
        end
    end
end
elem=elem(1:no,:);

%% FE-analysis
function [U]=FE(Nx,Ny,node,elem,x,penal)
[KE]=stiffness(elem);       % element stiffness matrix
Nn=size(node,1);    % Total number of nodes
Ne=size(elem,1);    % Total number of elements
K=zeros(3*Nn,3*Nn);     % global stiffness matrix
F=zeros(3*Nn,1);        % load vector
U=zeros(3*Nn,1);        % displacement vector
for jj=1:Ne
    n1=elem(jj,1);
    n2=elem(jj,2);
    edof=[3*n1-2;3*n1-1;3*n1;3*n2-2;3*n2-1;3*n2];
    K(edof,edof)=K(edof,edof)+x(jj)^penal*KE(:,:,jj);
end
% Define Loads and Supports
F(2,1)=-10;    %kN
fixeddofs=union([1:3*Nx:3*Nx*(Ny-1)+1],[3:3*Nx:3*Nx*(Ny-1)+3]);
fixeddofs=union([3*Nx-1],fixeddofs);
alldofs=[1:3*Nn];
freedofs=setdiff(alldofs,fixeddofs);
% Solving
U(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:);
U(fixeddofs,:)=0;

%% Stiffness matrix
function [KE]=stiffness(elem)
E=10000;
b=0.2;
h=0.4;
inertia=b*h^3/12;
csarea=b*h;
Ne=size(elem,1);    % Total number of elements
KE1=zeros(6,6,Ne);
KE2=zeros(6,6,Ne);
KE=zeros(6,6,Ne);
KEL=zeros(6,6,Ne);
for jj=1:Ne
    EL=elem(jj,3);
    KE1(:,:,jj)=inertia*E/(EL^3)*...
        [0  0       0       0   0       0;
        0   12      6*EL    0   -12     6*EL;
        0   6*EL    4*EL^2  0   -6*EL   2*EL^2;
        0   0       0       0   0       0;
        0   -12     -6*EL   0   12      -6*EL;
        0   6*EL    2*EL^2  0   -6*EL   4*EL^2];
    KE2(:,:,jj)=csarea*E/EL*...
        [1  0   0   -1  0   0;
        0   0   0   0   0   0;
        0   0   0   0   0   0;
        -1  0   0   1   0   0;
        0   0   0   0   0   0;
        0   0   0   0   0   0];
    cs=cos(elem(jj,4));
    si=sin(elem(jj,4));
    T=[cs si  0   0   0   0;
        -si cs  0   0   0   0;
        0   0   1   0   0   0;
        0   0   0   cs  si  0;
        0   0   0   -si cs  0;
        0   0   0   0   0   1];
    KEL(:,:,jj)=KE1(:,:,jj)+KE2(:,:,jj);
    KE(:,:,jj)=T'*KEL(:,:,jj)*T;
end

%% Optimality Criteria Update
function [xnew]=OC(x,elem,Vlim,dc)
l1=0; l2=100000; move=0.1;
while (l2-l1>1e-4)
    lmid=0.5*(l2+l1);
    xnew=max(0.001,max(x-move,min(1,min(x+move,x.*sqrt(-dc./lmid)))));
    vol=xnew'*elem(:,3);
    if vol-Vlim>0
        l1=lmid;
    else
        l2=lmid;
    end
end

Answer 1

TL;DR: answering the question directly

what this x.*sqrt(-dc./lmid) means, and why do I have to compare it?

(1) x.*sqrt(-dc./lmid) is x·B_e^η, the updated value x_e^new.

(2) You have to compare it to make sure that it is inside an appropriate interval (greater than or equal to x_min and x_e-m, and less than or equal to 1 and x_e+m).

Explanation

What is B_e^η

Following the optimality criteria, you want to update your previous x_e estimation to a more accurate value. You will achieve this by multiplying x_e by a factor q:

xenew = xe·q

• If q is exactly 1, then x_e is already in a (possibly local) optimum.
• If q > 1, then you underestimated x_e.
• If q < 1, then you overestimated x_e.

In your example, q is B_e^η.

Meaning of η

The first thing to note is that η is a tuning parameter, i.e., it is used to help the convergence of the method. It is called "damping coefficient" because it restricts the size of B_e^η (note that 0 < η <= 1). In your example, it was empirically set to 0.5, which is the square root.

Why the comparison

x_e is restricted to x_min <= x_e <= 1. So, we have to make sure that if x_e^new is less than x_min or greater than 1, we will truncate it. For reference, this operation is commonly called clamping. Besides, we also want to restrict x_e^new to be somewhat close to x_e (to help convergence). Thus, we need to clamp x_e^new again to be closer than an m value. Both comparisons will assure that x_e^new is valid (x_min <= x_e·B_e^η <= 1) and close to the previous value (x_e - m <= x_e·B_e^η <= x_e + m).

What is sqrt(-dc./lmid)

It is B_e^η. Looking at the image you added, we can see that the sqrt comes from η, -dc is -∂c/∂x_e, and lmid comes from λ·∂V/∂x_e. As in this case the volume of an element is always the same, dV = 1, and thus it gets out of the equation. Thus, lmid is only λ. It is called "mid" because to find it, the code runs a bisection algorithm, which divides the interval by half every iteration.