标准oc算法的推导与代码详解

对于变密度的参数化方法，设计变量x为材料相对密度，在已知材料的物性，包括弹性模型、密度以及给定载荷的条件下，我们希望简支梁的柔度最小，或者说使得结构势能最小（结构在力作用下的位移，该力做的功也就是势能）。那么当材料总体积保持为常数不变，给定载荷不变的情况下柔度最小可以理解为结构刚度最强。

问题描述

在结构力学相关设计中，通常会出现结构冗余的情况，例如受到一定载荷作用下的悬臂梁和简支梁其结构质量可以进一步减小。这一点在桥梁等设计中体现的更为明显。但是传统的桥梁设计方案是根据经验选择合适的桁架结构然后不断迭代设计方案。这种做法十分繁琐，目前我们可以通过拓扑优化^[1]的技术在概念设计阶段就提出最优方案，减少材料的浪费。

目前拓扑优化的参数化方法种类有很多，例如变密度法、水平集法、构型理论等。其中比较常见的是变密度法。密度法假设有一种材料密度可以在0-1之间变化。以这种参数为设计变量，密度为1代表材料存在，密度为0代表材料不存在。这样就在拓扑优化问题与材料分布问题之间建立了联系。

拓扑优化问题与传统最优化问题在基本原理上互通，其数学本质为求多元函数的极值问题。其基本形式可以进行如下抽象^[2]：

设$\mathbf{x}={\left(x_1,x_2,\cdot \cdot \cdot x_n\right)}^T$为n维欧氏空间${\mathbf{R}}^n$内一点，$f\left(\mathbf{x}\right),g_i\left(\mathbf{x}\right)\left(i=1,2,\cdot \cdot \cdot ,m\right)$和$h_j\left(\mathbf{x}\right)\left(j=m+1,\cdot \cdot \cdot ,p\right)$均为定的n元函数。最优化问题可表述为：在如下的约束条件下：
$$g_i\left(\mathbf{x}\right)\le 0i=1,2\cdot \cdot \cdot m$$

$$h_j\left(\mathbf{x}\right)=0j=m+1,\cdot \cdot \cdot p$$
求$f(x)$的最小值（或最大值）。通常称$f(x)$为目标函数，$g_i\left(\mathbf{x}\right)\le 0$为不等式约束条件，$h_j\left(\mathbf{x}\right)=0$为等式约束条件，$x$为设计变量。通用的表达形式为:

$$\{\begin{array}{c}\min f\left(\mathbf{x}\right)\\ s.t.g_i(\mathbf{x})\le 0i=1,2,\cdot \cdot \cdot m\\ h_j(\mathbf{x})=0j=m+1,\cdot \cdot \cdot p\\ x\in R^n\end{array}$$
上式又称优化列式，是待求解问题的表达式。
令:
R = { x | g i ( x ) ≤ 0 , i = 1 , 2 , ⋅ ⋅ ⋅ m ; h j ( x ) = 0 , j = m + 1 , ⋅ ⋅ ⋅ p } R = \left\{ x \middle| g_{i}\left( x \right) \leq 0,i = 1,2, \cdot \cdot \cdot m;h_{j}\left( x \right) = 0,j = m + 1, \cdot \cdot \cdot p \right\} R={x∣gi​(x)≤0,i=1,2,⋅⋅⋅m;hj​(x)=0,j=m+1,⋅⋅⋅p}
称$x\in R$
为上述问题的一个可行解集。

OC算法的数学描述

该问题可以写成如下数学形式：
$$\{\begin{array}{c}\underset{x}{\min }c\left(x\right)=U^{T}KU=\sum _{i=1}^N(x_i{)}^p{u}_i^T{k}_0{u}_i\\ s.t.\frac{V\left(x\right)}{V_0}=f\\ KU=F\\ 0<x_{min}\le x\le 1\end{array}$$
由于材料插值的缘故一般取材料的泊松比ν=1/3,相应的惩罚因子p=3。单元的刚度矩阵$k_i=x_i^k{k}_0$,其中$k_0$为单元刚度矩阵除以其弹性模量得到的“单位”刚度矩阵,结构的总刚矩阵为$K=\sum _{i=1}^N{k}_i$。有限元求解控制方程$U=K^{-1}F$可以得到每个结点的位移向量。在计算得到位移矩阵之后就可以通过标准优化准则（Standard Optimality Criteria，SOC)^[3]进行迭代。表达式如下：
$$x_i^{k+1}=\{\begin{array}{c}max(x_{min},x_i-m)if{x}_i{B}_i^{\eta }\le \left(x_{min,}{x}_i-m\right)\\ x_i{B}_i^{\eta },if\max (x_{min},x_i-m)<x_i{B}_i^{\eta }<min\left(x_{max},x_i+m\right)\\ min(x_{max,}{x}_i+m),ifmin\left(x_{max},x_i+m\right)\le x_i{B}_i^{\eta }\end{array}$$
其中$B_i=-\frac{\partial C}{\partial x_i}/\left(\frac{\partial V}{\partial x_i}\right)={p\left(x_i\right)}^{p-1}{u}_e^T{k}_0{u}_i/\left(v_e\right)$。

上式中m为一个移动极限因子,是为了保证算法的稳定性而引入的,在一次迭代中,一般不允许设计变量发生较大的变化;η(0<η<1)是阻尼系数,是为了使算法稳定,一般取0.5;$B_i$中的拉格朗日乘子l在每步迭代时,可以应用二分法根据结构的体积约束等式来确定。

当目标函数值的变化量足够小时停止迭代。

$$\left|\frac{C^{k+1}-C^k}{C^k}\right|<\epsilon 或\left|\frac{x^{k+1}-x^k}{x^k}\right|<\epsilon$$

OC法是拓扑优化常用方法，其通过定义Lagrange函数，将约束问题转化为非约束问题进行求解。

考虑约束条件$\frac{V(x)}{V_0}=f$，有Lagrange函数，形为：
$$L(x,\lambda )=c(x)+\lambda (V(x)-f{V}_0)$$
Karush-Kuhn-Tucker条件下，Lagrange函数的解满足：
$$\{\begin{array}{l}\frac{\partial L}{\partial x}=\frac{\partial c}{\partial x}+\frac{\lambda \partial V(r)}{\partial x}=0\\ \frac{\partial L}{\partial \lambda }=V(x)-f{V}_0=0\end{array}$$
Bendsøe提出的OC法过程包含两层循环，内层循环更新$x_e$以服从$\lambda$，外层循环更新$\lambda$以服从体积分数的约束。具体来说，OC法将在上述方程组第一等式成立条件下，定义以下尺度因子$B_i$：
$$B_i=-\frac{\frac{\partial c(x)}{\partial x_e}}{\frac{\lambda \partial V(x)}{\partial x_e}}$$
该尺度因子在$B_i=1$时满足第一等式。OC法根据$B_i$改变优化对象（$x_e$）的取值，方式如下：
$$x_e^{new}=x_e^{old}\sqrt{D_e},x_e^{min}\le x_e^{new}\le x_e^{max}$$
$B_i<1$意味着增大体积对减小柔度的效果比增大$x_e$取值的方式更强，应当减小$x_e$；反之，$D_e>1$时应当增大$x_e$。因为通常无法在单次迭代中得解，对$x_e$的限制是必要的；设计改变量$\Delta x_e$需满足$[x_e-∆x_{max},x_e+∆x_{max}]\subset [x_{min},x_{max}]$,使整体在设计范围内移动。

外层循环的更新方式各有不同，在本文介绍的拓扑优化方法中，作者使用二分法寻找Lagrange乘子$\lambda$。从二分集合[0，100000]开始，逐次减半范围，并且根据$\frac{\partial L}{\partial \lambda }=V(x)-f{V}_0$的符号，选取上下半集合作为新的二分集合，直到达到收敛条件。

结果展示

作者使用该代码得到如下结果：

OC算法的matlab代码及注释

该算法来自sigmond的99行代码^[3]
代码已上传至https://download.csdn.net/download/qq_42183549/74752801

%%%%%%%%  A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000  %%%%%%%%
%%%%%%%%  COMMENTED - OUT 1.0 BY HAOTIAN_W,                       JULY 2020  %%%%%%%%
function top(nelx,nely,volfrac,penal,rmin)
% ===================================================================================
% nelx    : 水平方向上的离散单元数;
% nely    : 竖直方向上的离散单元数;
%
% volfrac : 容积率，材料体积与设计域体积之比，对应的工程问题就是"将结构减重到百分之多少";
%
% penal   : 惩罚因子，SIMP方法是在0-1离散模型中引入连续变量x、系数p及中间密度单元，从而将离
%           散型优化问题转换成连续型优化问题，并且令0≤x≤1，p为惩罚因子，通过设定p>1对中间密
%           度单元进行有限度的惩罚，尽量减少中间密度单元数目，使单元密度尽可能趋于0或1;
% 
%           合理选择惩罚因子的取值，可以消除多孔材料，从而得到理想的拓扑优化结果:
%               X在迭代中最大置1，一旦置1即退出惩罚;
%               当penal<=2时     存在大量多孔材料，即可能存在停止迭代后依旧位于（0，1）间，散乱的元素，计算结果没有可制造性;
%               当penal>=3.5时   最终拓扑结果没有大的改变;
%               当penal>=4时     结构总体柔度的变化非常缓慢，迭代步数增加，计算时间延长，这是因为惩罚过重，惩罚加权后的元素敏度下降，
%               参见论文公式（4），dc的梯度包含Xe^penal的一项;
%            
% rmin    : 敏度过滤半径，防止出现棋盘格现象;棋盘格现象可看作是许多小尺度传力路径堆叠，敏度过滤限制传力路径的最小尺寸;
% 基于SIMP理论的优化准则法迭代分析流程:
%   1) 定义设计域，选择合适的设计变量、目标函数以及约束函数等其他边界条件;
%   2) 结构离散为有限元网格，计算优化前的单元刚度矩阵;
%   3) 初始化单元设计变量，即给定设计域内的每个单元一个初始单元相对密度;
%   4) 计算各离散单元的材料特性参数，计算单元刚度矩阵，组装结构总刚度矩阵，计算结点位移;
%   5) 计算总体结构的柔度值及其敏度值，求解拉格朗日乘子;
%   6) 用优化准则方法进行设计变量更新;
%   7) 检查结果的收敛性，如未收敛则转4)，循环迭代，如收敛则转8);
%   8) 输出目标函数值及设计变量值，结束计算。
% ===================================================================================
 
% x是设计变量（单元密度），给设计域内的每个单元一个初始相对密度，值为volfrac，即平均地减小各单元密度
x(1:nely,1:nelx) = volfrac; 
 
% loop储存迭代次数
loop = 0; 
 
% change储存每次迭代之后目标函数的改变值，用以判断是否收敛
change = 1.;
 
% 当目标函数改变量<=0.01时说明收敛，结束迭代
while change > 0.01
    loop = loop + 1; 
    
    % 将前一次的设计变量赋值给xold
    xold = x;  
 
    % 每次迭代都进行一次有限元分析，计算结点位移，并储存在全局位移数组U中，再从位移矩阵U中提取局部矩阵计算局部柔度
    [U] = FE(nelx,nely,x,penal);   
    
    % 计算单元刚度矩阵
    [KE] = lk;
    
    % c是用来储存目标函数的变量
    c = 0.;
    
    % 遍历设计域矩阵元素，从左上角第一个元素开始，一列一列
    for ely = 1:nely
        for elx = 1:nelx
            
            % 节点位移储存在U中，如果想获得节点位移必须先知道节点编号，进行索引;
            % 节点编号可以根据当前单元在设计域中的位置算出;
            % 每个元素均拥有四个节点，但仅需两个节点编号即可得知所有节点在位移矩阵中的位置，结合位移矩阵U作为单列向量的编号方式，需知左上与右上两个节点;
            % n1是左上角节点编号，n2是右上角节点编号;
            n1 = (nely+1)*(elx-1)+ely; 
            n2 = (nely+1)*elx+ely;
           
            % 局部位移数组Ue储存4个节点共8个自由度位移，每个节点分别有x、y两个自由度;，这是为什么需要为n1,n2乘以系数2
            % 因为是矩形单元，所以根据n1、n2两个节点的编号可以推演出单元所有节点的自由度编号;
            % 顺序是：[左上x；左上y；右上x；右上y；右下x；右下y；左下x；左下y];
            % 只适用于矩形单元划分网格;
            Ue = U([2*n1-1;2*n1; 2*n2-1;2*n2; 2*n2+1;2*n2+2; 2*n1+1;2*n1+2],1);
         
            % SIMP模型，将设计变量x从离散型变成指函数型，指数就是惩罚因子：x(ely,elx)^penal
            % 计算总体结构的柔度值，这里目标函数是柔度最小，参见论文中公式（1）
            c = c + x(ely,elx)^penal*Ue'*KE*Ue;
            
            % 计算总体结构的敏度值，实际上dc就是c对Xe的梯度，参见论文中公式（4）
            dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue;
        end
    end
    
    % 无关网格敏度过滤
    [dc] = check(nelx,nely,rmin,x,dc); 
    
    % 采用优化准则法（OC）求解当前模型，得出满足体积约束的结果，更新设计变量
    [x] = OC(nelx,nely,x,volfrac,dc);
    
    % 更新目标函数改变值，方式是各元素密度该变量的最大值小于目标（0.01）;
    change = max(max(abs(x-xold)));
    
    % 打印迭代信息: It.迭代次数，Obj.目标函数，Vol.材料体积比，ch.迭代改变量
    disp([' It.: ' sprintf('%4i',loop) ' Obj.: ' sprintf('%10.4f',c) ...
       ' Vol.: ' sprintf('%6.3f',sum(sum(x))/(nelx*nely)) ...
        ' ch.: ' sprintf('%6.3f',change )])
    
    % 当前迭代结果图形显示
    colormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6);
end
 
 
% 采用优化准则法（OC）迭代子程序 -------------------------------------------------------
% 数学模型主要求解算法有：优化准则法(OC)、序列线性规划法(SLP)、移动渐进线法(MMA); 
% OC适用于单约束优化问题求解，当求解复杂的多约束拓扑
% 优化问题时，采用SLP和MMA通常更方便;
% 参见论文中公式（2）（3）
function [xnew] = OC(nelx,nely,x,volfrac,dc)
% Input：  水平单元数nelx, 竖直单元数nely, 设计变量x, 材料体积比volfrac, 目标函数灵敏度dc;
% Output： 更新后的设计变量xnew;
 
% 定义一个取值区间，二分法，得到满足体积约束的拉格朗日算子
l1 = 0; l2 = 100000;
 
% 正向最大位移
move = 0.2;
 
while (l2-l1 > 1e-4)
    
    % 二分法，取区间中点
    lmid = 0.5*(l2+l1);
  
    % 参见论文中的公式（2）
    % sqrt(-dc./lmid)对应公式中Be^eta（eta=1/2），eta阻尼系数是为了确保计算的收敛性
    xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));
    
    % sum(sum(xnew))是更新后的材料体积, volfrac*nelx*nely是优化目标，用它们的差值判断是否收敛
    % sum()可以去看看help，注意一下行、列问题，不看也行，反正知道sum(sum())是所有元素求和就行
    if sum(sum(xnew)) - volfrac*nelx*nely > 0
        l1 = lmid;
    else
        l2 = lmid;
    end
end
 
 
% 无关网格敏度过滤子程序 --------------------------------------------------------------
% 参见论文中公式（5）（6）
function [dcn] = check(nelx,nely,rmin,x,dc)
% Input：  水平单元数nelx, 竖直单元数nely, 敏度过滤半径rmin, 设计变量x, 总体结构敏度dc;
% Output： 过滤后的目标函数敏度dcn;
 
% dcn清零，用来保存更新的目标函数灵敏度
dcn = zeros(nely,nelx);
 
for i = 1:nelx   
    for j = 1:nely
        sum=0.0; 
        
        % 在过滤半径定义的范围内遍历，遍历不超过元素边界
        for k = max(i-floor(rmin),1):min(i+floor(rmin),nelx)
            for l = max(j-floor(rmin),1):min(j+floor(rmin),nely)
                
                % 参见论文中公式（6）,fac即公式中卷积算子Hf
                % qrt((i-k)^2+(j-l)^2) 是计算此单元与相邻单元的距离，即公式中dist(e,f)，该卷积本质是加权平均
                %目的是将一定区域（rmin）内的元素联合计算敏度，避免小尺度传力路径出现
                fac = rmin-sqrt((i-k)^2+(j-l)^2);
                sum = sum+max(0,fac);
                
                % 参见论文中公式（5）
                dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k);
            end
        end
        dcn(j,i) = dcn(j,i)/(x(j,i)*sum);
    end
end
 
 
% 有限元求解子程序 --------------------------------------------------------------------
function [U] = FE(nelx,nely,x,penal)
% Input：  水平单元数nelx, 竖直单元数nely, 设计变量x, 惩罚因子penal;
% Output： 全局节点位移U;
 
% 计算单元刚度矩阵
[KE] = lk; 
 
% 总体刚度矩阵的稀疏矩阵
K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));
 
% 力矩阵的稀疏矩阵
F = sparse(2*(nely+1)*(nelx+1),1); 
 
% U清零，用来保存更新的全局节点位移
U = zeros(2*(nely+1)*(nelx+1),1);
 
for elx = 1:nelx
    for ely = 1:nely
        
        % 计算单元左上角、右上角节点编号
        n1 = (nely+1)*(elx-1)+ely; 
        n2 = (nely+1)* elx   +ely;
        
        % 同上主程序，计算单元4个节点8个自由度
        edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2];
        
        % 将单元刚度矩阵KE 组装成 总体刚度矩阵K
        K(edof,edof) = K(edof,edof) + x(ely,elx)^penal*KE;
    end
end
 
% 施加载荷，本算例应用了一个在左上角的垂直单元力
F(2,1) = -1;
 
% 施加约束，边界位移固定为0
% 消除线性方程中的固定自由度来实现支承结构，本算例左边第一列和右下角固定
fixeddofs = union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]);
 
% 剩下的不加约束的节点自由度，setdiff()从..中除去..
alldofs   = [1:2*(nely+1)*(nelx+1)];
freedofs  = setdiff(alldofs,fixeddofs);
 
% 求解线性方程组，得到各节点自由度的位移值储存在U中
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:); 
 
% 受约束节点固定自由度位移值为0
U(fixeddofs,:)= 0;
 
 
% 求解单元刚度矩阵子程序 --------------------------------------------------------------
% 有限元方法计算的一个重要的系数矩阵，表征单元体的受力与变形关系;
% 特点：对称性、奇异性、主对角元素恒正、奇数（偶数）行和为0;
% 矩形单元4节点 8*8矩阵;
function [KE] = lk
 
% 材料杨氏弹性模量
E = 1.; 
 
% 材料泊松比
nu = 0.3;
 
k=[ 1/2-nu/6   1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ... 
   -1/4+nu/12 -1/8-nu/8  nu/6       1/8-3*nu/8];
KE = E/(1-nu^2)*[ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)
                  k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)
                  k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)
                  k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)
                  k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)
                  k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)
                  k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)
                  k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];

参考文献

[1] Bendsøe Martin Philip, Kikuchi Noboru. Generating optimal topologies in structural design using a homogenization method[J]. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2):197-224.

[2] Bendsøe M，Sigmund O.Topology optimization theory，methods，and applications［M］.Germany: Springer Verlag，2003．

[3] O. Sigmund. A 99 line topology optimization code written in Matlab[J]. Structural and Multidisciplinary Optimization, 2001, 21(2) : 120-127.

ethod[J]. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2):197-224.

[4] 任晓辉. 连续体结构拓扑优化方法研究[D].长安大学,2007.