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