The Corresponding Files (Click to Save):

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:

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回圈，前兩層i1、j1是遍歷所有單元，后兩層i2、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行比大家都是構造矩陣，效率上沒有什么太大的差距，

% 單獨把這一段拎出來跑，測驗效率
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是所有單元第一個自由度，第一個自由度是啥，是左下角節點水平自由度，

接下來細說從edofVec到edofMat怎么回事，這里用到了repmat命令，具體去help，簡單點理解成復制就好了，

這里只要理解透 “一個節點所有自由度 = 該節點第一個自由度 + 相對值” 就結束了，

% 將列向量edofVec復制成8列，相當于所有單元8個自由度都初始化成各自的第一個自由度
repmat(edofVec,1,8)

% 在第一個自由度的基礎上加或減可以得到所有自由度
% 而且所有單元的操作是一樣的，因為是相對自身加減，單元大小、形狀相同，只有位置不同
repmat([0 1 2*nely+[2 3 0 1] -2 -1],nelx*nely,1)

對應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

轉載請註明出處，本文鏈接：https://www.uj5u.com/qianduan/52172.html

標籤：其他

上一篇：CSP202006-2 稀疏向量【map】

下一篇：資深大廠JAVA架構師帶你剖析Condition原始碼