The Corresponding Files (Click to Save):node

Code top88.m

References Efficient topology optimization in MATLAB using 88 lines of code

The Related Articals (Click in):算法

[TopOpt] 99行拓扑优化程序彻底注释

> [TopOpt] 针对99行改进的88行拓扑优化程序彻底注释

[TopOpt] 88行拓展的82行拓扑优化程序彻底注释. PDE

[TopOpt] 88行拓展的71行拓扑优化程序彻底注释. CONV2

The program can be promoted by line:spring

top88(120,40,0.5,3.0,3.5,1)

代码注释

88行程序为了提升效率，逻辑、流程没有99行程序那么清楚，建议初学先读99行。数组

%%%%%%%%%%%% AN 88 LINE TOPOLOGY OPTIMIZATION CODE              Nov 2010 %%%%%%%%%%%%
%%%%%%%%%%%% COMMENTED - OUT BY HAOTIAN_W                    AUGUST 2020 %%%%%%%%%%%%
function top88(nelx,nely,volfrac,penal,rmin,ft)
% ===================================================================================
% 88行程序在99行程序上的主要改进:
%   1) 将for循环语句向量化处理，发挥MATLAB矩阵运算优点;
%   2) 为不断增长数据的数组预分配内存，避免MATLAB花费额外时间寻找更大的连续内存块;
%   3) 尽量将部分程序从循环体里抽出，避免重复计算;
%   4) 设计变量再也不表明单元伪密度，新引入真实密度变量xphys;
%   5) 将原先的全部子程序都集成在主程序里，避免频繁调用;
%   整体上，程序的效率有显著提高（近百倍）、内存占用下降，可是对初学者来讲可读性不如99行
% ===================================================================================
% nelx    : 水平方向上的离散单元数;
% nely    : 竖直方向上的离散单元数;
%
% volfrac : 容积率，材料体积与设计域体积之比，对应的工程问题就是"将结构减重到百分之多少";
%
% penal   : 惩罚因子，SIMP方法是在0-1离散模型中引入连续变量x、系数p及中间密度单元，从而将离
%           散型优化问题转换成连续型优化问题，而且令0≤x≤1，p为惩罚因子，经过设定p>1对中间密
%           度单元进行有限度的惩罚，尽可能减小中间密度单元数目，使单元密度尽量趋于0或1;
% 
%           合理选择惩罚因子的取值，能够消除多孔材料，从而获得理想的拓扑优化结果:
%               当penal<=2时     存在大量多孔材料，计算结果没有可制造性;
%               当penal>=3.5时   最终拓扑结果没有大的改变;
%               当penal>=4时     结构整体柔度的变化很是缓慢，迭代步数增长，计算时间延长;
%            
% rmin    : 敏度过滤半径，防止出现棋盘格现象;

% ft      : 与99行程序不一样的是，本程序提供了两种滤波器，
%           ft=1时进行灵敏度滤波，获得的结果与99行程序同样；ft=2时进行密度滤波。
% ===================================================================================
%% 定义材料属性
% E0杨氏弹性模量;
E0 = 1;

% Emin自定义空区域杨氏弹性模量，为了防止出现奇异矩阵，这里和99行程序不一样，参见论文公式（1）;
% 不要定义成0，不然就没有意义了;
Emin = 1e-9;

% nu泊松比;
nu = 0.3;

%% 有限元预处理
% 构造平面四节点矩形单元的单元刚度矩阵KE，详见有限元理论推导
A11 = [12  3 -6 -3;  3 12  3  0; -6  3 12 -3; -3  0 -3 12];
A12 = [-6 -3  0  3; -3 -6 -3 -6;  0 -3 -6  3;  3 -6  3 -6];
B11 = [-4  3 -2  9;  3 -4 -9  4; -2 -9 -4 -3;  9  4 -3 -4];
B12 = [ 2 -3  4 -9; -3  2  9 -2;  4  9  2  3; -9 -2  3  2];
KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]);

% nodenrs存放单元节点编号，按照列优先的顺序，从1到(1+nelx)*(1+nely);
nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx);

% edofVec存放全部单元的第一个自由度编号（左下角），参见下面图1;
edofVec = reshape(2*nodenrs(1:end-1,1:end-1)+1,nelx*nely,1);

% edofMat按照行存放每一个单元4个节点8个自由度编号，因此列数是8，行数等于单元个数;
% 存放顺序是：[左下x  左下y  右下x  右下y  右上x  右上y  左上x  左上y];
% 第一个repmat将列向量edofVec复制成8列，其余自由度从第一个自由度上加或减能够获得，参见下面图2;
edofMat = repmat(edofVec,1,8)+repmat([0 1 2*nely+[2 3 0 1] -2 -1],nelx*nely,1);

% 根据iK、jK和sK三元组生成整体刚度矩阵的稀疏矩阵K，K = sparse(iK,jK,sK)
iK = reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1);
jK = reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1);

% 施加载荷，直接构造稀疏矩阵
F = sparse(2,1,-1,2*(nely+1)*(nelx+1),1);

% 施加约束，与99行程序相同，惟一的区别是这里把“定义载荷约束”放在了循环体外，提升效率
U = zeros(2*(nely+1)*(nelx+1),1);
fixeddofs = union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]);
alldofs = [1:2*(nely+1)*(nelx+1)];
freedofs = setdiff(alldofs,fixeddofs);

%% 滤波预处理
% 根据iH、jH和sH三元组生成加权系数矩阵的稀疏矩阵H，H = sparse(iH,jH,sH);
% H是论文公式（8）中的Hei,Hs是论文公式（9）中的Sigma(Hei);

% 为数组iH、jH、sH预分配内存;
iH = ones(nelx*nely*(2*(ceil(rmin)-1)+1)^2,1);
jH = ones(size(iH));
sH = zeros(size(iH));
k = 0;

% 4层for循环，前两层i一、j1是遍历全部单元，后两层i二、j2是遍历当前单元附近的单元
% 这种敏度过滤技术的本质是利用过滤半径范围内各单元敏度的加权平均值代替中心单元的敏度值
for i1 = 1:nelx
  for j1 = 1:nely
    e1 = (i1-1)*nely+j1;
    for i2 = max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx)
      for j2 = max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely)
        e2 = (i2-1)*nely+j2;
        k = k+1;
        iH(k) = e1;
        jH(k) = e2;
        sH(k) = max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2));
      end
    end
  end
end
H = sparse(iH,jH,sH);
Hs = sum(H,2);

%% 迭代初始化
x = repmat(volfrac,nely,nelx);   % x设计变量;
xPhys = x;                       % xphys单元物理密度，真实密度，这里与99行不同;
loop = 0;                        % loop存放迭代次数;
change = 1; 

%% 进入优化迭代，到此为止上面的部分都是在循环外，比99行效率提升不少
while change > 0.01
  loop = loop + 1;
  
  %% 有限元分析求解
  % (Emin+xPhys(:)'.^penal*(E0-Emin))就是论文公式（1），由单元密度决定杨氏弹性模量;
  sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(E0-Emin)),64*nelx*nely,1);
  
  % 组装整体刚度矩阵的稀疏矩阵;
  % K = (K+K')/2确保整体刚度矩阵是彻底对称阵，由于这会影响到MATLAB求解有限元方程的算法;
  % 当K是实对称正定矩阵时则采用Cholesky平方根分解法，反之则采用速度更慢的LU三角分解法;
  K = sparse(iK,jK,sK); K = (K+K')/2;
  
  % 正式求解,KU=F;
  U(freedofs) = K(freedofs,freedofs)\F(freedofs);
  
  %% 目标函数和关于真实密度场的灵敏度信息
  % 参见论文公式（2），ce是ue^T*k0*ue，c是目标函数;
  ce = reshape(sum((U(edofMat)*KE).*U(edofMat),2),nely,nelx);
  c = sum(sum((Emin+xPhys.^penal*(E0-Emin)).*ce));
  
  % 参见论文公式（5），只不过将设计变量x换成真实密度xphys;
  dc = -penal*(E0-Emin)*xPhys.^(penal-1).*ce;
  
  % 参见论文公式（6）;
  dv = ones(nely,nelx);
  
  %% 敏度滤波 或 密度滤波
  
  if ft == 1
    % ft=1灵敏度滤波（获得的结果同99行程序），参见论文公式（7）;
    % 完成敏度滤波;
    % 1e-3是公式（7）中的gamma，max(1e-3,x(:))是为了防止分母出现0;
    dc(:) = H*(x(:).*dc(:))./Hs./max(1e-3,x(:));
    
  elseif ft == 2
    % ft=2密度滤波
    % 密度滤波进行两个操做，一个是密度滤波，在下面的循环里面
    % 另外一个是根据链式法则修正目标函数和体积约束的灵敏度信息，就是这里，参见公式（10）;
    dc(:) = H*(dc(:)./Hs);
    dv(:) = H*(dv(:)./Hs);
  end
  
  %% OC优化准则法更新设计变量和单元密度
  l1 = 0; l2 = 1e9; move = 0.2;
  while (l2-l1)/(l1+l2) > 1e-3
    lmid = 0.5*(l2+l1);
    
    % 更新设计变量，参见论文公式（3）;
    xnew = max(0,max(x-move,min(1,min(x+move,x.*sqrt(-dc./dv/lmid)))));
    
    if ft == 1
      % 敏度滤波没有密度滤波那么复杂，设计变量就是当前单元的伪密度;
      xPhys = xnew;
      
    elseif ft == 2
      % 完成密度滤波，参见论文公式（9），设计变量通过滤波以后才是单元伪密度;
      xPhys(:) = (H*xnew(:))./Hs;
    end
    if sum(xPhys(:)) > volfrac*nelx*nely, l1 = lmid; else l2 = lmid; end
  end
  change = max(abs(xnew(:)-x(:)));
  x = xnew;
  
  %% 显示结果（同99行程序）
  fprintf(' It.:%5i Obj.:%11.4f Vol.:%7.3f ch.:%7.3f\n',loop,c, ...
    mean(xPhys(:)),change);
  colormap(gray); imagesc(1-xPhys); caxis([0 1]); axis equal; axis off; drawnow;
end

论文公式汇总

Modified SIMP方法中基于单元伪密度的杨氏弹性模量	$E_{e}\left(x_{e}\right)=E_{\min }+x_{e}^{p}\left(E_{0}-E_{\min }\right), \quad x_{e} \in[0,1]$	（1）
优化问题（柔度最小化）的数学表述	$\begin{array}{cl} \min _{\mathbf{x}}: & c(\mathbf{x})=\mathbf{U}^{\mathrm{T}} \mathbf{K} \mathbf{U}=\sum_{e=1}^{N} E_{e}\left(x_{e}\right) \mathbf{u}_{e}^{\mathrm{T}} \mathbf{k}_{0} \mathbf{u}_{e} \\ \text { subject to: } & V(\mathbf{x}) / V_{0}=f \\ & \mathbf{K U}=\mathbf{F} \\ & 0 \leq \mathbf{x} \leq \mathbf{1} \end{array}$	（2）
OC优化准则法更新设计变量	$x_{e}^{\mathrm{new}}=\left\{\begin{array}{ll} \max \left(0, x_{e}-m\right) & \text { if } x_{e} B_{e}^{\eta} \leq \max \left(0, x_{e}-m\right) \\ \min \left(1, x_{e}+m\right) & \text { if } x_{e} B_{e}^{\eta} \geq \min \left(1, x_{e}-m\right) \\ x_{e} B_{e}^{\eta} & \text { otherwise } \end{array}\right.$	（3）
	$B_{e}=\frac{-\frac{\partial c}{\partial x_{e}}}{\lambda \frac{\partial V}{\partial x_{e}}}$	（4）
目标函数关于设计变量的灵敏度信息	$\frac{\partial c}{\partial x_{e}}=-p x_{e}^{p-1}\left(E_{0}-E_{\min }\right) \mathbf{u}_{e}^{\mathrm{T}} \mathbf{k}_{0} \mathbf{u}$	（5）
体积约束关于设计变量的灵敏度信息	$\frac{\partial V}{\partial x_{e}}=1$	（6）
对目标函数的灵敏度信息进行敏度滤波	$\frac{\widehat{\partial c}}{\partial x_{e}}=\frac{1}{\max \left(\gamma, x_{e}\right) \sum_{i \in N_{e}} H_{e i}} \sum_{i \in N_{e}} H_{e i} x_{i} \frac{\partial c}{\partial x_{i}}$	（7）
权重系数矩阵	$H_{e i}=\max \left(0, r_{\min }-\Delta(e, i)\right)$	（8）
对设计变量进行密度滤波获得单元伪密度	$\tilde{x}_{e}=\frac{1}{\sum_{i \in N_{e}} H_{e i}} \sum_{i \in N_{e}} H_{e i} x_{i}$	（9）
在密度滤波中，根据链式法则修正目标函数和体积约束关于设计变量的灵敏度信息	$\frac{\partial \psi}{\partial x_{j}}=\sum_{e \in N_{j}} \frac{\partial \psi}{\partial \tilde{x}_{e}} \frac{\partial \tilde{x}_{e}}{\partial x_{j}}=\sum_{e \in N_{j}} \frac{1}{\sum_{i \in N_{e}} H_{e i}} H_{j e} \frac{\partial \psi}{\partial \tilde{x}_{e}}$	（10）

四节点矩形单元刚度矩阵

% 单元刚度矩阵
A11 = [12  3 -6 -3;  3 12  3  0; -6  3 12 -3; -3  0 -3 12];
A12 = [-6 -3  0  3; -3 -6 -3 -6;  0 -3 -6  3;  3 -6  3 -6];
B11 = [-4  3 -2  9;  3 -4 -9  4; -2 -9 -4 -3;  9  4 -3 -4];
B12 = [ 2 -3  4 -9; -3  2  9 -2;  4  9  2  3; -9 -2  3  2];
KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]);

这里把有限元单元刚度矩阵的基础知识本身推导一遍就基本没问题了，网上不少资料，这里很少说了。这一段跟99行比你们都是构造矩阵，效率上没有什么太大的差距。app

% 单独把这一段拎出来跑，测试效率
99行
Elapsed time is 0.000402 seconds.
88行
Elapsed time is 0.000366 seconds.

单元节点自由度编号

% 分别是节点编号、单元第一个自由度编号、全部自由度编号
nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx);
edofVec = reshape(2*nodenrs(1:end-1,1:end-1)+1,nelx*nely,1);
edofMat = repmat(edofVec,1,8)+repmat([0 1 2*nely+[2 3 0 1] -2 -1],nelx*nely,1);

这里就拿4X3网格举例子了。程序注释上面有，很绕口，不如直观点来看看这三个矩阵到底在干吗。注意nodenrs是节点编号，不是单元编号，nelx个单元一条边有(nelx+1)个节点这没什么好解释的。edofVec是全部单元第一个自由度，第一个自由度是啥，是左下角节点水平自由度。ide

图1 矩形单元4节点8自由度编号矩阵（矢量化）

接下来细说从edofVec到edofMat怎么回事。这里用到了repmat命令，具体去help，简单点理解成复制就行了。函数

这里只要理解透 “一个节点全部自由度 = 该节点第一个自由度 + 相对值” 就结束了。oop

% 将列向量edofVec复制成8列，至关于全部单元8个自由度都初始化成各自的第一个自由度
repmat(edofVec,1,8)

% 在第一个自由度的基础上加或减能够获得全部自由度
% 并且全部单元的操做是同样的，由于是相对自身加减，单元大小、形状相同，只有位置不一样
repmat([0 1 2*nely+[2 3 0 1] -2 -1],nelx*nely,1)

图2 edofMat的矢量化构造思路

对应99行程序中 Line: 17-21遍历单元，当时的思路是遍历每一个单元计算一下左上角、右上角节点编号，而后推算该单元全部自由度编号。测试

而在这里，88行程序中，做者使用矢量化思路，经过矩阵运算来提高程序运行效率。从上面咱们也看到了，因为全部的单元大小形状都是同样的只有位置不同，只要定了每一个单元第一个自由度编号（图2 等号左边第一个矩阵），剩余的操做都是同样的（图2 等号左边第二个矩阵）。优化

下面是我写的一段测试代码，分别把99行和88行中相应片断粘贴过来，稍微修改使两段程序获得相同的结果。而后测试处理4万网格二者花费的时间。

从结果看，处理4万网格，矢量化程序运行时间不到遍历单元方法的万分之二，差别十分显著。何况，在99行程序中，每一个迭代循环都要跑一次这个，若是有300次优化迭代那就要跑300次；而在88行程序中总共只须要在开头进行一次。这也就是为何我用3700X处理器跑3万网格的99行程序要花掉我大半个下午，而一样3万网格的88行程序喝口水就能跑完。并且屡次测试验证，网格越多，效率差距越大。

% 网格编号处理部分运行效率测试
clear
clc
nelx = 2e2;
nely = 2e2;
% 99行程序 遍历单元
disp('遍历单元')
tic,
edof1 = [];
for elx = 1:nelx
    for ely = 1:nely
        n1 = (nely+1)*(elx-1)+ely; 
        n2 = (nely+1)*elx+ely;
        edof1 = [edof1;[2*n1+1 2*n1+2 2*n2+1 2*n2+2 2*n2-1 2*n2 2*n1-1 2*n1]];
    end
end
toc
% 88行程序 矢量化
disp('矢量化')
tic,
nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx);
edofVec = reshape(2*nodenrs(1:end-1,1:end-1)+1,nelx*nely,1);
edof2 = repmat(edofVec,1,8)+repmat([0 1 2*nely+[2 3 0 1] -2 -1],nelx*nely,1);
toc

% i5 7400 @3.0GHz
遍历单元
Elapsed time is 21.128882 seconds.
矢量化
Elapsed time is 0.003363 seconds.

% Ryzen7 3700X @4.2GHz
遍历单元
Elapsed time is 12.857031 seconds.
矢量化
Elapsed time is 0.002607 seconds.

有限元求解

% 循环体外面
iK = reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1);
jK = reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1);

% 循环体里面
sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(E0-Emin)),64*nelx*nely,1);
K = sparse(iK,jK,sK); K = (K+K')/2;
U(freedofs) = K(freedofs,freedofs)\F(freedofs);

根据iK、jK和sK三元组生成整体刚度矩阵的稀疏矩阵K，索引向量iK和jK已经在循环体外面定义过了，sK须要在循环内肯定。sK根据单元刚度矩阵KE和单元杨氏弹性模量求得，单元杨氏弹性模量参见论文公式（1）。

组装整体刚度矩阵的稀疏矩阵。K = (K+K')/2 是为了确保整体刚度矩阵是彻底对称阵，由于这会影响到MATLAB求解有限元方程时使用的算法，当K是实对称正定矩阵时MATLAB采用“Cholesky平方根分解法”求解，若是K不是对称正定则采用速度更慢的“LU三角分解法”。

下面这段测试程序很好地展现了LU分解和对称正定矩阵的Cholesky分解在MATLAB中运行效率的差距。不难发现，88行程序始终围绕着“高效”，做者也为此下了不少功夫。

% 比较MATLAB中LU分解和Cholesky分解的运行速度
clear,clc
A = gallery('lehmer',1e4);
disp('LU分解')
tic, lu(A); toc
disp('Cholesky分解')
tic, chol(A); toc

LU分解
Elapsed time is 5.887353 seconds.
Cholesky分解
Elapsed time is 2.854781 seconds.

密度滤波和敏度滤波

if ft == 1
    % 完成敏度滤波
    dc(:) = H*(x(:).*dc(:))./Hs./max(1e-3,x(:));
  elseif ft == 2
    % 修正灵敏度值
    dc(:) = H*(dc(:)./Hs);
    dv(:) = H*(dv(:)./Hs);
  end

xnew = max(0,max(x-move,min(1,min(x+move,x.*sqrt(-dc./dv/lmid)))));
    if ft == 1
      xPhys = xnew;
    elseif ft == 2
      % 这里才真正完成密度滤波
      xPhys(:) = (H*xnew(:))./Hs;
    end

88行程序提供了两种滤波方法，指定ft=1时使用敏度滤波，指定ft=2时使用密度滤波。

注意这里其实分了两段，这是由于密度滤波并非一次性完成的，在循环体外先根据链式法则对目标函数和体积约束的灵敏度信息进行修正，而后循环体内对设计变量进行密度滤波获得单元伪密度。

而敏度滤波就要简单的多，设计变量就表明单元伪密度，直接对目标函数的灵敏度信息进行滤波就能够。

在密度滤波中，引入了“设计变量（Design Variables）”和“物理密度（Physical Density）”两个场，在程序中的符号分别是x和xphys。这种状况下，设计变量并不表明真实的单元伪密度场，通过滤波的物理密度才是真实的单元伪密度场。包括以后的OC优化准则法迭代过程当中，拉格朗日算子lmid也是由真实物理密度xphys决定的。在推导目标函数/约束条件关于设计变量的灵敏度信息时，必须使用链式法则来肯定，参见论文公式（10）。

使用88行MATLAB拓扑优化程序获得的MBB梁优化结果对比
（a）60 X 20；（b）150 X 50；（c）300 X 100
上面三幅图像是使用敏度滤波获得的结果，下面三幅图像是使用密度滤波获得的结果

程序运行测试

读完了99行程序和针对99行改进的88行程序，就让咱们一块儿来对比测试一下吧。把下面这段程序加进主程序先后，正常运行一下就获得结果啦。我用一个90*30的网格就在笔记本上简单跑了一遍，结果来看88行程序在效率上的提高仍是很显著的，几乎只用了99行六分之一的时间。

tic

% 主程序段

t = toc;
Mem = memory;
Mem.MemUsedMATLAB = Mem.MemUsedMATLAB/1e6;
disp(['Iterations    :   ' sprintf('%4i',loop) ' times'])
disp(['Elapsed Time  :    ' sprintf('%6f',t) ' s'])
disp(['Memory Used   :    ' sprintf('%6.2f',Mem.MemUsedMATLAB) ' MB'])

% 99行程序
Iterations     :    379 times
Objective Func :    196.5764
Elapsed Time   :    64.179272 s
Memory Used    :    1683.04 MB

% 88行程序
Iterations     :    221 times
Objective Func :    197.6266
Elapsed Time   :    12.723110 s
Memory Used    :    1690.55 MB

虽然在前面咱们已经提到过好几回这段程序提升MATLAB运行效率的方法，这里咱们再总结一下吧：

尽可能减小for循环，培养矢量化思路；
为数组预分配内存空间；
尽可能选用合适的MATLAB函数，这一点很宽泛，须要不少积累才行，这里咱们至少知道了一条，Cholesky比LU快；
尽可能避免频繁调用子程序；
精简循环体，能放在外面的代码就放到循环外面。

关于MATLAB的效率优化我一直想写一篇博客，但是我懒，涉及到太多方面一直不想动手写，hhh我尽快。

参考资料

[1] Andreassen, E., Clausen, A., Schevenels, M. et al. Efficient topology optimization in MATLAB using 88 lines of code. Struct Multidisc Optim 43, 1–16 (2011).

本文连接：http://blog.csdn.net/BAR_WORKSHOP/article/details/108287668