-
TopOpt | 99行拓扑优化程序完全注释
2020-08-31 00:47:13:)Click Here For More Details. 最近学习Sigmund的99行Matlab拓扑优化程序,参考了其他大神们的注释,作为初学者还有不少看不懂的地方,结合论文学习又补充整理了一遍。博客搬家到自己搭建的 主页 啦q(≧▽≦q),大家快来逛逛鸭!
The Corresponding Files (Click to Save):
- Code: top.m
- References: A 99 line topology optimization code written in MATLAB
The Related Artical (Click in):
The program can be promoted by line:
top(30,10,0.5,3.0,1.5)
代码注释
以前学了点皮毛就直接用商业软件搬砖了,总觉得有点心虚,所以入门第一步,拜读了Sigmund的99行Matlab拓扑优化程序,参考了其他大神们的注释,作为初学者还有不少看不懂的地方,我自己又补充整理了一遍。我是零基础小白,所以代码前面先交代了一些理论背景,我自己能看懂相信所有人都能看懂,哈哈哈。
%%%%%%%% 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; % % 合理选择惩罚因子的取值,可以消除多孔材料,从而得到理想的拓扑优化结果: % 当penal<=2时 存在大量多孔材料,计算结果没有可制造性; % 当penal>=3.5时 最终拓扑结果没有大的改变; % 当penal>=4时 结构总体柔度的变化非常缓慢,迭代步数增加,计算时间延长; % % rmin : 敏度过滤半径,防止出现棋盘格现象; % =================================================================================== % 结构优化的数学模型常用名词: % 1) 设计变量 在设计中可调整的、变化的基本参数,本算例是单元密度; % 2) 目标函数 设计变量的函数,优化设计的目标,本算例是柔度最小; % 3) 约束条件 几何、应力、位移等,难点在建立约束方程,本算例是几何体积约束; % 4) 终止准则 结束迭代的条件,本算例是目标函数变化量<=0.010; % 5) 载荷工况 定义结构所有可能的受力情况,本算例是单一载荷; % =================================================================================== % 基于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,x用来储存这一次的结果,之后还要比较它们以判断是否收敛 xold = x; % 每次迭代都进行一次有限元分析,计算结点位移,并储存在全局位移数组U中 [U] = FE(nelx,nely,x,penal); % 计算单元刚度矩阵 [KE] = lk; % c是用来储存目标函数的变量 c = 0.; % 遍历设计域矩阵元素,从左上角第一个元素开始,一列一列 for ely = 1:nely for elx = 1:nelx % 节点位移储存在U中,如果想获得节点位移必须先知道节点编号,进行索引; % 节点编号可以根据当前单元在设计域中的位置算出; % n1是左上角节点编号,n2是右上角节点编号; n1 = (nely+1)*(elx-1)+ely; n2 = (nely+1)*elx+ely; % 局部位移数组Ue储存4个节点共8个自由度位移,每个节点分别有x、y两个自由度; % 因为是矩形单元,所以根据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); % 更新目标函数改变值 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) 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; % 施加约束,消除线性方程中的固定自由度来实现支承结构,本算例左边第一列和右下角固定 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) OC准则优化更新 (2) (3) 目标函数灵敏度 (4) 灵敏度滤波 (5) 卷积算子(加权因子) (6) 单元刚度矩阵
程序里最后一部分子程序,详细请学习有限元的基础知识:平面四节点矩形单元的单元刚度矩阵推导。
参考资料
版权声明:本文为博主原创文章,转载请附上原文出处链接和本声明。
本文链接:http://blog.csdn.net/BAR_WORKSHOP/article/details/108274360
-
电机优化程序使用说明
2017-02-08 09:32:43电机优化程序使用说明电机优化程序使用说明 引言 1编写目的 2项目背景 3 定义 软件概述 1目标 2功能 3 性能 运行环境 1硬件 硬件最小配置 硬件推荐配置 2支持软件 使用说明 1安装和初始化 第一步解压缩软件包到一个...电机优化程序使用说明 1.引言
1.1编写目的
本手册主要针对电机设计人员而编写,在文中主要介绍了电机优化程序的安装、初始化和使用方法。对于不是很熟悉电机设计和ANSYS Maxwell使用的人员并不适用。
1.2项目背景
电机优化程序是在香港理工大学电磁实验室傅为农教授的指导下,由多名研究人员一起合作完成的项目。其中,吴会欢负责算法实现、优化流程、程序界面设计和文档编写。盛田田、翁旭、林启芳、马悦等人负责测试以及反馈错误。
1.3 定义
在电机优化程序中,定义了如下默认设置:最优方向为目标函数趋向于0
2. 软件概述
2.1目标
希望可以通过一个程序,实现对电动机结构的快速优化,提高电动机的性能,不仅好,而且快。(Computer Automated Design)
2.2功能
电机优化程序是一个使用多目标遗传算法优化电机的程序,它可以配合ANSYS Maxwell优化电机设计。
2.3 性能
a.数据精确度
所有参数最高保留6位小数。
b.时间特性
对于每个任务的设置阶段,程序相应应该在1秒左右,一般不会假死。如果遇到假死问题请关闭程序并重新打开。
c.灵活性
电机优化程序使用ANSYS Maxwell作为评估性能的程序,如果ANSYS Maxwell没有安装在C盘默认路径,请修改PATH。
例如:版本:16.1
C:\Program Files\AnsysEM\AnsysEM16.1\Win64
上述路径需要添加在PATH中,如果将ANSYS Maxwell安装在E盘的a文件夹中(
E:\a
),那么需要添加的路径则是:E:\a\AnsysEM\AnsysEM16.1\Win64
3. 运行环境
3.1硬件
硬件最小配置:
CPU: 1核
内存大小: 1G
硬盘容量: 20G硬件推荐配置:
CPU: 16核或更多
内存大小: 16GB或更多
硬盘容量: 100GB或更多3.2支持软件
操作系统:Windows 7或者更高版本, 64bit
支持软件:ANSYS Maxwell 16.02或ANSYS Maxwell 16.1, 64bit
注意:本软件不支持ANSYS Electromagnet Suite 17
4. 使用说明
4.1安装和初始化
第一步:解压缩软件包到一个目录下。
直接使用压缩软件解压即可。
第二步:设置系统PATH变量
将maxwell.exe所在的目录加入PATH
对于Maxwell安装在程序默认路径的用户来说,可以参照一下设置:版本:16.1
C:\Program Files\AnsysEM\AnsysEM16.1\Win64
用户可以将上述灰色背景的路径复制到资源管理器,然后查看是否存在maxwell.exe这个文件。如果存在的话将上述路径添加至PATH中。
注意:在将路径复制到PATH中时不可以删除原有字符,需要在保留之前字符串,然后在之后输入一个英文分号
;
后加入C:\Program Files\AnsysEM\AnsysEM16.1\Win64
这段字符串第三步:打开gaWizard.exe
双击gaWizard.exe即可
4.2输入
注意:单目标优化暂时不支持,请使用Maxwell自带的遗传算法实现单目标优化。
电机优化程序中,用户需要准备一个ANSYS Maxwell设计文件,文件格式为mxwl
,例如motor.mxwl
。
目前,电机优化程序可以支持电机的转矩(Torque)、效率(Efficiency)、功率因数(Power factor)、反电动势(Back EMF)、转矩波动(Ripple)优化,这些参数可以是2个或者多个。在本说明中,我们需要利用这些参数设计目标函数来实现多目标优化。其他的一些参数列在下面:
名称 属性 支持的电机类型 无限制 目标函数数量 无限制,一般不超过5个 变量数量 无限制,一般不超过50个 4.2.1数据背景
用户需要准备的数据包括电机设计中需要改变的变量名称,变量的范围(上下限)以及单位。
用户可以准备一个下图这样的单子:4.2.2数据格式
在mxwl文件中,请准备一个design且仅有一个。如下图:
对于优化的目标函数,请在Output Variables里设置,如下图
4.2.3输入举例
在文件目录下,有一个叫
2.mxwl
的例子文件。用户可以打开它研究一下需要设置的地方应该如何设置。几个关键的部分
1. Design Properities
2. Output Variables
例如在Output Variables中,torque ripple表示为ripple,它使用的是Moving1.Torque的数据,然后使用ripple函数获得ripple。
3. Ripple图表的设置
上图中ripple的Y轴数据来自ripple(Moving1.Torque),而不是直接的转剧torque。
4. Power Factor图表的设置
power factor的计算需要电流和反电动势,所以在新建报告的时候需要选择电流和电压两部分。如下图所示:
设置完成便如下图所示:4.3输出
电机优化程序会输出3个文件,OBJ存放最终的目标函数的结果,VAR存放变量和目标函数的结果,用|符号分隔开。4.3.1数据背景
OBJ、VAR和result.txt文件将会放在程序所在目录下。4.3.2数据格式
4.3.3举例
4.4出错和恢复
如何判断出错
- 程序在运行,但是打开windows任务管理器后发现CPU使用率长时间小于5%,且maxwell.exe进程远小于CPU核心数。例如8核电脑运行程序一段时间后,进度条一直没有前进,而且打开任务管理器后发现只剩下1个或者2个maxwell.exe进程,且相关maxwellEngine和solver2d进程也无反应。
- 解决方法:
这个问题通常是由于ANSYS Maxwell 程序假死造成,可以通过打开任务管理器,关闭maxwell.exe进程。(一般只有1个或者2个maxwell.exe进程)
- 解决方法:
- 在有进度条的界面中,进度条飞快结束
- 解决方法:
PATH设置有误,请重新设置PATH环境变量。
- 解决方法:
4.5求助查询
如果遇到无法解决问题请将出错截图和具体的问题描述发送至whhxp1028@gmail.com 吴会欢
5. 运行说明
5.1运行表
5.2运行步骤
5.2.1运行控制
5.2.2操作信息
5.2.3输入/输出文件
5.2.4启动或恢复过程
6. 非常规过程
7. 操作命令一览表
8. 程序文件(或命令文件)和数据文件一览表
名称 文件名 主程序名 gaWizard.exe 优化程序名 jmetaltest.exe 结果保存文件名 result.txt 默认配置文件名 generalSettings.ini 9. 用户操作举例
第一步:准备Maxwell模型,以永磁同步电机为例子
在优化程序的目录中有一个new.mxwl作为例子。
首先用户需要画一个永磁同步电机。
例如下图的电机:然后在Design Properities设置变量
然后设置目标函数,用于导出
右键点击Results,选择Output Variables打开后如下图所示
这里取”efficiency”, “powerfactor”, “ripple”, “torque” 和 “voltage”作为优化的目标函数。由于遗传算法定义目标函数的值趋向0时为最优,所以一些越大越好的参数可以通过一个参考值减去实际值来调整,例如可以用1-Pout/Pin来实现efficiency的目标函数。完成设置之后便如上图。
接下来需要创建Report用于导出。请注意,对于每个目标函数都要建立一个Transient Report。
选择Results点击右键,选择Create Transient Report中的Rectangular Report,分别建立Report。新建完成后请将Report的名字命名为相应的目标函数的名字(注意大小写)。建立完成后如下图所示
注意:如果需要优化powerfactor请在powerfactor表中包含2个值,分别是同相的Induced Voltage和InputCurrent。需要优化ripple时要使用range function中的ripple函数处理torque。
第二步:设置环境变量
设置PATH
优化程序需要调用maxwell 进行求解,所以用户需要设置windows环境变量。
首先找到控制面板,打开高级系统设置,如下图所示:点击右下方环境变量按钮,出现如下对话框:
在系统变量中找到PATH这个值,如果没有可以自己新建一下。
点击编辑,出现下图对话框,不要将对话框中选中的字符串删除
在字符串的最后加入ANSYS Maxwell的路径,例如在本机上,ANSYS Maxwell 16.1中的maxwell.exe安装在
C:\Program Files\AnsysEM\AnsysEM16.1\Win64
所以将上述字符串复制,并且粘贴到最后。如果之前字符串最后没有分号,请手工输入英文分好将两个字符串隔开。完成后如下图所示:
点击确定关闭对话框。然后点击剩下的对话框的确定按钮关闭所有对话框,完成环境变量配置。
第三步:使用优化程序
新用户取得程序时,是下图中的一个压缩包:
解压缩
解压缩后,会有一个相应文件名的文件夹出现,里面有一个release文件夹
打开release文件夹,里面有下图中的文件:
其中,gaWizard.exe就是主程序,双击打开即可。
打开后,画面如图所示:首先需要激活程序,根据对话框提示,点击help按钮,这时会跳出一个对话框,上面显示如下信息:
上图中有一个机器码,它已经被复制到剪切板中。直接点击OK按钮关闭对话框。
这时候,请将机器码粘贴,发送至相关人员处,然后获取类似格式的激活码。
接下来,点击对话框下列中最左边的Activate按钮,输入激活码。如果错误,则会提示
那么请重新找相关人员进行反馈。
如果正常激活,则会提示然后点击Next进入优化程序主页面。如果正常激活,那么之后打开程序都将直接进入主页面,如下图
主页面中,用户可以选择2种方式调用优化程序,使用一步一步的方式设置配置文件,或者使用现有的配置文件直接进行优化。
对于新用户,请选择第一种,第一个配置文件设置界面如下图图所示。点击File Path 列最右边的
...
按钮,打开之前准备的电机模型。打开后,程序将自动读取design的名称,以及之前在Output Variables里面定义的目标函数的名称。
接下来选择需要优化的目标函数,这里选择torque,powerfactor和efficiency。
然后,点击下一步,选择变量,这里用户一共定义了16个参数,我们选择其中几个参数作为变量,并且设置上限、下限以及单位,STEP功能尚未完工。
完成设置后点击Next按钮,进入遗传算法设置页面。
如下图:这里基本不用更改,如果只是测试,请将offspring,population和maxEvaluations的值改小一些,例如offsprings和population改成10,maxEvaluations改成40.
然后点击Next,这里用户需要查看之前的设置是否有错,还可以保存或者打开配置文件。如下图所示。如果没错,请点击I Confirm…,然后点击Next按钮进入优化环节。
开始优化页面就很简单了,点击start开始任务
然后等待进度条跑完即可。结束之后,点击下一步进入查看结果页面。
结果页面显示的是刚刚结束计算的一代的结果。|
号前面显示的是变量,变量名称按照英文字母顺序排列|
号后面显示的是目标函数,目标函数的名称按照字母顺序排序 - 程序在运行,但是打开windows任务管理器后发现CPU使用率长时间小于5%,且maxwell.exe进程远小于CPU核心数。例如8核电脑运行程序一段时间后,进度条一直没有前进,而且打开任务管理器后发现只剩下1个或者2个maxwell.exe进程,且相关maxwellEngine和solver2d进程也无反应。
-
Sigmund的99行Matlab拓扑优化程序简析
2018-05-29 12:34:42引言 Sigmund在2001年在Structural and Multidisciplinary Optimization 上...该论文后附带了一个Matlab拓扑优化程序。这个只有99行代码程序基于Matlab环境构建了一个完整的拓扑优化流程,其中包含:前处理(构建...引言
Sigmund在2001年在Structural and Multidisciplinary Optimization 发表一篇名为 “A 99 line topology optimization code written in Matlab”论文。该论文附带一个使用Matlab编写的拓扑优化程序。这个只有99行代码的程序基于Matlab环境构建了一个完整的拓扑优化流程:前处理(构建有限元仿真模型), 有限元模型分析计算,拓扑优化迭代和后处理(分析结果显示)。Sigmund在论文中对这个程序的拓扑优化流程做了详细解释。本文仅从程序设计角度解析这段代码。
Sigmund的99行Matlab拓扑优化程序使用模块化方法设计,主要包含以下几个模块:
- 程序主流程
- 有限元模型求解模块
- Filter模块
- 单元刚度阵计算模块:计算平面四边形单元的刚度矩阵
- 优化模块: 使用优化准则法更新设计变量
- 目标函数计算和灵敏度分析模块
下面几个小节首先列出Sigmund99行拓扑优化程序完整代码,接着解释程序中5个全局变量在优化流程中的功用,然后对程序的各个功能模块做详细解析,最后总结程序的局限性并列出扩展程序功能时所需要修改的相应模块。99行拓扑优化代码
Sigmund的99行Matlab拓扑优化程序如下所示
% a 99 line topology optimization code by Ole Sigmund,October 1999 clear nelx=60; nely=40; volfrac=0.5; penal=3.; rmin=1.5; % initialize x(1:nely,1:nelx)=volfrac; loop=0; change=1; % start ineration while change>0.01 loop=loop+1; xold=x; % FE analysis [U]=FE(nelx,nely,x,penal); % objective function and sensitivity analysis [KE]=lk;; c=0.; for ely=1:nely for elx=1:nelx n1=(nely+1)*(elx-1)+ely; n2=(nely+1)*elx +ely; 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); c=c+x(ely,elx)^penal*Ue'*KE*Ue; dc(ely,elx)=-penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue; end end % filtering of sensitivities [dc]=check(nelx,nely,rmin,x,dc); % design update by the optimality criteria method [x]=oc(nelx,nely,x,volfrac,dc); % print result change=max(max(x-xold)) disp(['It.:' sprintf( '%4i',loop) ' Obj.:' sprintf(' %10.4f',c) ... ' Vol.:' sprintf('%6.3f',sum(sum(x))/(nelx*nely)) ... ' ch.:' sprintf('%6.3f',change)]) % plot densities colormap(gray);imagesc(-x);axis equal;axis tight; axis off;pause(1e-6); end % FE analysis function [U]=FE(nelx,nely,x,penal) [KE]=lk; K=sparse(2*(nelx+1)*(nely+1),2*(nelx+1)*(nely+1)); F=sparse(2*(nely+1)*(nelx+1),1); U=sparse(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; edof=[2*n1-1;2*n1;2*n2-1;2*n2;2*n2+1;2*n2+2;2*n1+1;2*n1+2]; K(edof,edof)=K(edof,edof)+x(ely,elx)^penal*KE; end end % define loads and supports ip=(nelx+1)*(nely+1); F(2*ip,1)=-1; fixeddofs =[1:2*(nely+1)]; alldofs =[1:2*(nely+1)*(nelx+1)]; freedofs =setdiff(alldofs,fixeddofs); % solving U(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:); U(fixeddofs,:)=0; % mesh-independency filter function [dcn]=check(nelx,nely,rmin,x,dc) 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) fac=rmin-sqrt((i-k)^2+(j-l)^2); sum=sum+max(0,fac); 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 % Element stiffness matrix 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)]; % optimality criteria update function [xnew]=oc(nelx,nely,x,volfrac,dc) l1=0; l2=100000; move=0.2; 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))))); if sum(sum(xnew))-volfrac*nelx*nely>0; l1=lmid; else l2=lmid; end end
变量初始化
nelx=60; % x轴方向单元数量 nely=40; % y轴方向单元数量 volfrac=0.5; % 单元的材料体积比 penal=3.; % 惩罚因子,将设计变量从[0,1]转化为指数函数, % 主要目的是使单元材料密度更加"黑白分明", % 惩罚因子就是这个指数函数的指数, 其值越大,效果越明显 rmin=1.5; % 过滤半径,目的是预防出现棋盘格现象
设计变量初始化
x(1:nely,1:nelx)=volfrac; % 所有设计变量的赋初始值
主程序
主程序迭代流程
% start ineration while change>0.01 loop=loop+1; % loop是迭代计数 xold=x; % 将设计变量保存在 变量xold中 % FE analysis [U]=FE(nelx,nely,x,penal); % 有限元模型分析,将计算得到的各个节点位移值保存在数组U中 % objective function and sensitivity analysis [KE]=lk; % 计算单元刚度阵 c=0.; % 保存目标函数(柔度)的变量, % 遍历所有单元 for ely=1:nely for elx=1:nelx % 根据单元在设计域中位置 计算单元节点编号 n1=(nely+1)*(elx-1)+ely; % 左上角节点编号 n2=(nely+1)*elx +ely; % 右上角节点编号 % 根据以上两个节点的编号 推演出 单元所有节点的自由度编号 % 这就限制了此程序只能使用四边形单元剖分四边设计域所得到的有限元模型 % % 从数组 U 中根据节点自由度编号 提取该单元的 位移向量 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(Solid Isotropic Material with Penalization)材料插值模型 c=c+x(ely,elx)^penal*Ue'*KE*Ue; % 计算灵敏度 dc(ely,elx)=-penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue; end end % filtering of sensitivities [dc]=check(nelx,nely,rmin,x,dc); % design update by the optimality criteria method [x]=oc(nelx,nely,x,volfrac,dc); % print result change=max(max(x-xold)) disp(['It.:' sprintf( '%4i',loop) ' Obj.:' sprintf(' %10.4f',c) ... ' Vol.:' sprintf('%6.3f',sum(sum(x))/(nelx*nely)) ... ' ch.:' sprintf('%6.3f',change)]) % plot densities colormap(gray);imagesc(-x);axis equal;axis tight; axis off;pause(1e-6); end
优化模型为:
代码 x(ely,elx)^penal*Ue’*KE*Ue 就是计算
对 的梯度是
此程序使用的SIMP(Solid Isotropic Material with Penalization)材料插值模型,具体代码为:
x(ely,elx)^penal ;网格划分为
1-------(nely+1)+1------------------------------------nely+1)*nelx | (1,1) | 2-------(nely+1)+2 | ... -----------(nely+1)*(elx-1)+ely----(nely+1)*elx+ely | | (elx,ely) | ... | nely+1---(nely+1)*2----------
单元节点编号顺序为:
Node_1 Node_2 (nely+1)*(elx-1)+ely----(nely+1)*elx+ely | (elx,ely) | (nely+1)*(elx-1)+ely+1----(nely+1)*elx+ely+1 Node_4 Node_3
过滤器模块
%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%% function [dcn]=check(nelx,nely,rmin,x,dc) % 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) % sqrt((i-k)^2+(j-l)^2) 是计算此单元与相邻单元的距离 fac = rmin-sqrt((i-k)^2+(j-l)^2); sum = sum + max(0,fac); 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
Sigmund的论文中对此函数有详细解释,论文中符号与代码中变量名的对应关系如下表所示:
变量 论文中符号 rmin fac dc dcn 优化迭代模块
% optimality criteria update function [xnew]=oc(nelx,nely,x,volfrac,dc) l1=0; l2=100000; move=0.2; 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))))); if sum(sum(xnew))-volfrac*nelx*nely>0; l1=lmid; else l2=lmid; end end
这段代码就是实现下面这个方程式:
有限元模型求解模块
% FE analysis function [U]=FE(nelx,nely,x,penal) [KE]=lk; %计算单元刚度矩阵 K=sparse(2*(nelx+1)*(nely+1),2*(nelx+1)*(nely+1)); F=sparse(2*(nely+1)*(nelx+1),1); U=sparse(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; edof=[2*n1-1;2*n1;2*n2-1;2*n2;2*n2+1;2*n2+2;2*n1+1;2*n1+2]; K(edof,edof)=K(edof,edof)+x(ely,elx)^penal*KE; end end % define loads and supports ip=(nelx+1)*(nely+1); F(2*ip,1)=-1; % 施加载荷 % 施加位移约束 fixeddofs =[1:2*(nely+1)]; alldofs =[1:2*(nely+1)*(nelx+1)]; % 所有不受约束的节点自由度 freedofs =setdiff(alldofs,fixeddofs); % solving 求解线性方程组, 得到节点自由度的位移值 U(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:); U(fixeddofs,:)=0; % 受约束的节点自由度的位移值 设为 0
单元刚度阵模块
% Element stiffness matrix 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)];
4节点矩形单元的4个节点在单元局部坐标系中坐标分别为:
Node x y 0 0 1 0 1 1 0 1 各节点形函数设为
以各节点形函数为基函数构建线性单元的位移场。单元内任一点 处的位移可由四个节点的变形通过单元形函数求得:
其中 分别为各个节点在单元局部坐标系中 坐标轴向和 坐标轴向的位移分量。
根据几何方程,应变为:
则单元B阵为
材料本构阵为:
单元刚度阵为:
其中
其余的依次类推。程序的局限性
Sigmund的99行Matlab拓扑优化程序 主要关注于拓扑优化算法是验证. 以牺牲了程序的通用性为代价,短短的99行代码中完成一个完整的拓扑优化流程。整个程序代码简洁、计算高效,但各程序模块的功能能省就省,往往使用最简单最直接的方法实现,存在一定的局限性,例如:
- 设计域默认是矩形
- 只能处理一种单元– 平面四边形单元
- 只能使用一种各向同性材料
- Filter函数依赖于设计域网格划分方式修改程序
如果要改变优化模型和有限元模型,调整项与相应的功能模块如下表所示:
调整项 功能模块 函数名 优化模型 主程序 优化算法 目标函数计算模块 OC 边界条件 有限元模型求解模块 FE 材料类型 单元刚度阵模块 lk 过滤算法 过滤器模块 check 材料插值模型 主程序 -
CSAPP:优化程序性能
2010-12-10 23:04:00优化程序性能这章作为CSAPP中最闪光的一章,其重要程度不言而喻。此实验分为了两个部分:第一部分是对一个多项式计算的优化;第二部分是对矩阵代码的优化。 首先,我们必须了解一些优化程序的常识。 编写高效的...优化程序性能这章作为CSAPP中最闪光的一章,其重要程度不言而喻。此实验分为了两个部分:第一部分是对一个多项式计算的优化;第二部分是对矩阵代码的优化。
首先,我们必须了解一些优化程序的常识。
编写高效的程序需要两个方面下足马力:第一,我们必须选择一组最好的的算法和数据结构;第二,我们必须编写出编译器能够有效优化以转换成高效可执行代码的源代码。那么此时,我们就需要理解优化编译器的能力和局限性。它的局限性表现在两个方面:第一,存储器别名;
提示:考虑*xp == *yp的情况。
第二,函数调用。
提示:考虑fun(x)里面是一个全局变量或静态变量,并且递增。
接下来将逐个介绍优化程序的方法。
第一个方法,代码移动。什么是代码移动呢?它就是把在循环中执行多次但结果却固定不变的代码移动到循环之外,这是最简单的优化。
第二个方法,减少过程调用。为什么要减少函数调用呢?因为每次调用一个函数,这个函数都可能会执行冗余的操作,比如数组边界检查,类型检查等。举例来说,如果我们需要调用一个函数来获得数组中的一个元素,我们可以首先获得数组首地址,然后通过索引访问每个元素,而非调用一个函数,然后检查数组是否越界,然后才返回元素值。
第三种方法,消除不必要的存储器引用。也就是说,在循环体中不使用指针来存放结果,而使用临时变量。为什么呢?目的就是减少读、写内存地址的时间。
随着对优化的深入,研究汇编代码以及理解现代处理器是如何执行计算则变得很重要。下面继续了解几种优化技术。
第四种方法,循环展开 。什么是循环展开?其思想是在一次迭代中对多个数组元素进行访问和合并操作。那这样做有什么原因或者说带来了什么好处呢?受限于资源约束,每个周期内只有两个功能单元可以执行数据相关的操作,然而一个周期内可能有涉及到数据的多个操作。那么,这些操作就不能在同一周期内完成,必须暂停等待,降低性能。通过循环展开可以在每次迭代中执行更多的数据操作,避免一个周期内产生多于两个操作,从而降低循环开销。理想情况下的整数加法调度
实际资源受限的情况下,整数加法的调度
优化后,三次循环展开的整数加法的调度
第五种方法,迭代分割 。什么是迭代分割呢?我们知道处理器的几个功能单元是流水线化的,这说明它们在一个操作完成之前可以开始另一个新的操作。但是,我们的代码没有这种能力。当执行当前操作时,处理器会处理其它事务(一般不会暂停),直到本操作完成,才开始下一个操作。这严重影响了程序的性能。而通过迭代分割,则可以很好地缓解此问题,以提高性能。迭代分割其实就是把一个操作分割成两个或更多的部分,并在最后合并结果以提高性能。产生二路并行或n路并行,提高了程序的并行性。
理想情况下的整数乘法调度
实际资源受限的情况下,整数乘法的调度
优化后,二次展开、二次并行的整数乘法操作调度
比较直观的例子
运用上面的方法来优化第一部分多项式代码,思考!
对于第二部分的矩阵优化,我们用到的核心思想就是利用局部性原理,包括时间局部性和空间局部性,避免大幅度地跨越内存块。具体的示例在第三章“程序的机器级表示”的嵌套数组以及第六章“存储器层次结构”的局部性原理中讲述。
-
oracle优化程序模式
2005-09-22 09:10:00优化程序模式• 基于规则的:– 使用 --- 套等级系统– 语法驱动和数据字典驱动• 基于成本的:– 选择成本最低的路径– 统计数据驱动基于规则的优化在基于规则的优化中通过检查查询服务器进程选择它访问数据的路径该... -
优化程序性能 计算机系统结构 深入理解计算机系统
2015-10-23 14:33:58优化程序性能 1编写高效程序:合适的数据结构和算法,编译器能够有效优化以转换为高效可执行代码的 源码,对处理量特别大的计算将任务分为多个部分; 程序优化:消除不必要的内容(函数调用,条件测试,存储器... -
优化程序性能总结
2015-10-23 14:34:57性能优化有三个层次: 系统层次 算法层次 代码层次 系统层次关注系统的控制流程和数据流程,优化主要考虑如何减少消息传递的个数;如何使系统的负载更加均衡;如何充分利用硬件的性能和设施;如何减少系统额外... -
《深入理解计算机系统》优化程序性能
2013-12-22 16:05:45优化程序性能的几点: 1)高级设计。选择恰当的数据结构和算法; 2)基本编码原则:消除连续的函数调用;消除不必要的存储器引用; 3)低级优化。循环展开;提高并行性,如多个累积变量,重新结合技术,条件传送代替... -
使用数据库连接池优化程序性能
2011-11-21 17:55:08使用数据库连接池优化程序性能 在我们做项目时连接数据库一般采用两种方法:1、是应用程序直接获取连接;2、是通过数据库连接池来回去连接。 第一种方法是用户每次请求都需要向数据库获得链接,这样有一个很大... -
VTune 分析和优化程序性能的工具
2009-02-24 15:53:00Intel VTune Performance Analyzer是一个用于分析和优化程序性能的工具,它能确定你的程序的hotspot,帮助你找到导致性能不理想的原因,从而让你能据此对程序进行优化。 1. What is VTune? Intel VTune Performance ... -
多目标优化程序测试程序(VC++)
2010-05-27 07:48:00多目标优化测试程序。 -
CSAPP:优化程序性能(一)
2017-06-20 13:39:56程序员必须在实现和维护程序的简单性和运算速度之间做出权衡,几分钟就能编写一个简单的插入程序,而一个高效的排序算法程序可能需要一天或更长时间来实现和优化, 大多数编译器,例如GCC向用户提供了一些对它们所... -
伪谱法优化程序已经开发出来
2017-08-09 22:29:10可以用于弹道优化等场合。《最优控制问题的Legendre...伪谱法程序编起来比我想象得要困难很多,编了整整两个月,克服了大大小小的困难。 目前优化算法采用snopt优化工具箱,我觉得遗传算法等也可以,后续将进行试验。 -
第五章 优化程序性能
2015-10-11 11:45:51写程序的最主要目标就是使它在所有可能的情况下都正确...2. 编写出编译器能够有效优化以转换成高效可执行代码的源代码。 3. 针对运算量特别大的计算,并行计算(将一个任务分成多个部分,这些部分可以在多核和多处 -
最优化程序c++的方法
2015-08-27 01:53:50解决的问题最优化问题的一般形式 minf(paras)min \quad f(paras)s.t.paras∈[paras_lower,paras_upper]s.t. \quad paras\in[paras\_lower,paras\_upper] 这个是数学上的一般形式,当求函数的最大值时候只要加上... -
性能之巅:定位和优化程序CPU、内存、IO瓶颈
2020-12-09 14:29:50不同程序有不同的性能关注点,比如科学计算关注运算速度,游戏引擎注重渲染效率,而服务程序追求吞吐能力。 服务器一般都是可水平扩展的分布式系统,系统处理能力取决于单机负载能力和水平扩展能力,所以,提升单机... -
【优化程序性能】编译器的优化和限制
2014-05-02 23:00:11写程序最主要的目标是使他在所有的可能的情况下dou -
深入理解计算机系统:优化程序性能
2016-09-16 22:34:37编写高效程序需要几类活动: 选择合适的算法和数据结构; 编写出编译器能够有效优化以转换成高效可执行代码的源代码; 针对运算量特别大的计算,将一个任务分成多个部分,在多核和多处理器的某种组合上并行地计算。 ... -
Matlab粒子群算法(PSO)优化程序——经典实例
2019-03-31 10:29:47粒子群算法( Particle Swarm Optimization, PSO)最早是由Eberhart和Kennedy于1995年提出,它的基本概念源于对鸟群觅食行为的研究...(1)粒子:优化问题的候选解,指鸟群中的一个个个体; (2)位置:候选解所在... -
MATLAB 优化程序【profile简明用法】
2018-05-09 22:45:49编程过程中,经常需要评估哪一部分代码比较耗费时间,这对于优化代码非常重要。Visual Studio中的profile功能可以评测,Matlab同样也有这个功能,而且使用起来也比较简单。-基本命令 profile on : 开启profile ... -
深入理解计算机系统之旅(五)优化程序性能
2014-05-28 11:42:36好的算法和数据结构在编写高性能的程序时固然重要,但是却不是全部,如果我们想要写出性能更好的程序就需要了解编译器是否如何工作和优化我们的代码的,当然并不是要求所有的程序员都去了解和掌握此技能,仅仅对程序... -
TopOpt | 针对99行改进的88行拓扑优化程序完全注释
2020-09-14 20:39:23The 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行拓扑优化程序完全注释 ...
-
个人发卡系统无后门.zip
-
浅析facebook的信息架构
-
JMETER 性能测试基础课程
-
用户管理-源码
-
用于非侵入式负载监控(能量分解)的迁移学习
-
角型打字稿-源码
-
Cookie&Session笔记
-
边缘计算白皮书.rar
-
【动态规划】Lintcode 114. 不同的路径
-
MMM 集群部署实现 MySQL 高可用和读写分离
-
SE-ML-Algorithms-DataStruc:在软件工程与机器学习中尝试一些简单的练习,算法,数据结构和软件模式。 语言会有所不同,包括Python,Java,CC ++,R,Julia,Golang,Haskell,Scala,Javascript,Matlab和SQL-源码
-
敏捷开发工具ScrumWorks使用简介
-
Galera 高可用 MySQL 集群(PXC v5.6 + Ngin
-
基于Qt的LibVLC开发教程
-
2021 年该学的 CSS 框架 Tailwind CSS 实战视频
-
Vue watch与watchEffect的区别
-
libFuzzer视频教程
-
PPT大神之路高清教程
-
使用vue搭建微信H5公众号项目
-
综合城市能源系统的最佳投资计划.zip