博客搬家到自己搭建的 主页 啦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):

>TopOpt | 99行拓扑优化程序完全注释
	
TopOpt | 针对99行改进的88行拓扑优化程序完全注释
	
TopOpt | 88行拓展的82行拓扑优化程序完全注释. PDE
	
TopOpt | 88行拓展的71行拓扑优化程序完全注释. CONV2


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）
		
单元刚度矩阵

程序里最后一部分子程序，详细请学习有限元的基础知识：平面四节点矩形单元的单元刚度矩阵推导。



参考资料

[1] Sigmund, O. A 99 line topology optimization code written in Matlab. Struct Multidisc Optim 21, 120–127 (2001). 

版权声明：本文为博主原创文章，转载请附上原文出处链接和本声明。

本文链接：http://blog.csdn.net/BAR_WORKSHOP/article/details/108274360

电机优化程序使用说明


电机优化程序使用说明
引言
1编写目的
2项目背景
3 定义

软件概述
1目标
2功能
3 性能

运行环境
1硬件
硬件最小配置
硬件推荐配置

2支持软件

使用说明
1安装和初始化
第一步解压缩软件包到一个目录下
第二步设置系统PATH变量
第三步打开gaWizardexe

2输入
21数据背景
22数据格式
23输入举例
几个关键的部分

3输出
31数据背景
32数据格式
33举例
4出错和恢复
5求助查询


运行说明
1运行表
2运行步骤
21运行控制
22操作信息
23输入输出文件
24启动或恢复过程


非常规过程
操作命令一览表
程序文件或命令文件和数据文件一览表
用户操作举例
第一步准备Maxwell模型以永磁同步电机为例子
第二步设置环境变量
第三步使用优化程序







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开始任务 


然后等待进度条跑完即可。

结束之后，点击下一步进入查看结果页面。 

结果页面显示的是刚刚结束计算的一代的结果。

|号前面显示的是变量，变量名称按照英文字母顺序排列 
|号后面显示的是目标函数，目标函数的名称按照字母顺序排序

引言

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 

优化模型为：



 $$\left\{
 \begin{array}{c c c l }
 min & c(x) & = & U^{T}KU = \sum_{e = 1}^{n} (x_{e})^{p} u_{e}^{T}K_{e} u_{e}   \\
 s.t. & & \\
  & \frac{V_{\rho}}{V} & = & f  \\
  & KU & = & F  \\
  & x_{min} & < & x < 1        \\
 \end{array}
 \right.$$ 


代码 x(ely,elx)^penal*Ue’*KE*Ue 就是计算 $(x_{e})^{p} u_{e}^{T}K_{e} u_{e}$ 

对$x_{e}$ 的梯度是 

 $$\frac{dc}{dx} = - p (x_{e})^{p-1} u_{e}^{T}K_{e} u_{e}$$  

此程序使用的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
  
$r_{min}$


  
fac
  
$H_{f}$


  
dc
  
$\frac{c}{x_{e}}$


  
dcn
  
$\bar {  \frac{c}{x_{e}} }$






优化迭代模块



% 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             

这段代码就是实现下面这个方程式: 

 $$x_{e}^{new} =
\left \{
\begin{array}{ c c c c l}
max(x_{min}, x_{e} - m) & if &                        &      & x_{e} B_{e}^{\eta} \leq max(x_{min}, x_{e}- m)  \\
x_{e} B_{e}^{\eta}      & if & max(x_{min}, x_{e}- m) & \leq & x_{e} B_{e}^{\eta} \leq min(1, x_{e}+m ) \\
min(1, x_{e}+m )        & if &   min(1, x_{e}+m )     & \leq & x_{e} B_{e}^{\eta} \\
\end{array}
\right  .$$ 




有限元模型求解模块



% 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



  
$(x_{1}, y_{1})$
  
0
  
0


  
$(x_{2}, y_{2})$
  
1
  
0


  
$(x_{3}, y_{3})$
  
1
  
1


  
$(x_{4}, y_{4})$
  
0
  
1




各节点形函数设为 

 $$\left \{
    \begin{array}{c c c c c}
 N_{1}  & =   &  (x - 1 ) (y - 1)  \\  
 N_{2}  & =   & - x (y - 1)  \\ 
 N_{3}  & =  &  x y  \\
 N_{4}  & =  &  -(x - 1 ) y  \\
    \end{array}
  \right .$$  

以各节点形函数为基函数构建线性单元的位移场。单元内任一点$(x,y)$ 处的位移可由四个节点的变形通过单元形函数求得： 

 $$\left \{
    \begin{array}{c c l}
u(x,y) & = & \sum_{i=1}^{4}N_{i}(x,y)u_{i} \\
v(x,y) & = & \sum_{i=1}^{4}N_{i}(x,y)v_{i} \\
    \end{array}
  \right  .$$  

其中$u_{i}, v_{i}$ 分别为各个节点在单元局部坐标系中$x$ 坐标轴向和$y$ 坐标轴向的位移分量。 

根据几何方程，应变为： 

 $$\xi_{xx} = \frac{\partial u }{\partial x} = \frac{\partial N_{1}}{\partial x} u_{1} + \cdots + \frac{\partial N_{4}}{\partial x} u_{4}\\
$$ 




 $$\xi_{yy} = \frac{\partial v}{\partial y} = \frac{\partial N_{1}}{\partial y} v_{1} + \cdots + \frac{\partial N_{4}}{\partial y} v_{4} \\
$$ 




 $$\xi_{xy} = \frac{\partial v}{\partial x} + \frac{\partial u }{\partial y} 
      = \frac{\partial N_{1}}{\partial y} u_{1} + \frac{\partial N_{1}}{\partial x} v_{1} + \cdots 
       + \frac{\partial N_{4}}{\partial y} u_{4} + \frac{\partial N_{4}}{\partial x} v_{4}$$   


则单元B阵为 

 $$B =
\left [
    \begin{array}{c c  c}
\frac{\partial N_{1}(x,y)}{\partial x} & \cdots & 0 \\
  0  & \cdots & \frac{\partial N_{4}(x,y)}{\partial y}   \\
\frac{\partial N_{1}(x,y)}{\partial y}   & \cdots & \frac{\partial N_{4}(x,y)}{\partial x}  \\
    \end{array}
  \right ]$$ 


材料本构阵为： 

 $$D
 =
\frac{E}{1-\nu^{2}}
\left [
    \begin{array}{c c c}
1  &  \nu  & 0 \\
\nu & 1 & 0  \\
0 &  0  & \frac{1-\nu}{2} \\
    \end{array}
  \right ]$$ 


单元刚度阵为： 

 $$K = \int B^{T} D B ds$$  

其中 

 $$K_{11} =\frac{E}{1-\nu^{2}} \int_{0}^{1}\int_{0}^{1} ( (y-1)^{2} + \frac{1-\nu}{2}(x-1)^{2} )dxdy 
= \frac{E}{1-\nu^{2}}( \frac{1}{2} - \frac{ \nu }{6})$$  

其余的$K_{ij}$依次类推。




程序的局限性

Sigmund的99行Matlab拓扑优化程序 主要关注于拓扑优化算法是验证.   以牺牲了程序的通用性为代价，短短的99行代码中完成一个完整的拓扑优化流程。整个程序代码简洁、计算高效，但各程序模块的功能能省就省，往往使用最简单最直接的方法实现，存在一定的局限性，例如： 

 - 设计域默认是矩形 

 - 只能处理一种单元– 平面四边形单元  

 - 只能使用一种各向同性材料 

 - Filter函数依赖于设计域网格划分方式



修改程序

如果要改变优化模型和有限元模型，调整项与相应的功能模块如下表所示：




  
调整项
  
功能模块
  
函数名



  
优化模型
  
主程序
  



  
优化算法
  
目标函数计算模块
  
OC


  
边界条件
  
有限元模型求解模块
  
FE


  
材料类型
  
单元刚度阵模块
  
lk


  
过滤算法
  
过滤器模块
  
check


  
材料插值模型
  
主程序
  

优化程序性能这章作为CSAPP中最闪光的一章，其重要程度不言而喻。此实验分为了两个部分：第一部分是对一个多项式计算的优化；第二部分是对矩阵代码的优化。


首先，我们必须了解一些优化程序的常识。

编写高效的程序需要两个方面下足马力：第一，我们必须选择一组最好的的算法和数据结构；第二，我们必须编写出编译器能够有效优化以转换成高效可执行代码的源代码。那么此时，我们就需要理解优化编译器的能力和局限性。它的局限性表现在两个方面：
第一，存储器别名；
void twiddle1(int *xp,int *yp)
{
	*xp += *yp;
	*xp += *yp;
}
void twiddle2(int *xp,int *yp)
{
	*xp += 2 * (*yp);
}  
提示：考虑*xp == *yp的情况。 
第二，函数调用。
int fun(int);
int func1(x)
{	
	fun(x)+fun(x)+fun(x)+fun(x);
}
int func2(x)
{
	4*fun(x);
}
提示：考虑fun(x)里面是一个全局变量或静态变量，并且递增。 

接下来将逐个介绍优化程序的方法。

第一个方法，代码移动。什么是代码移动呢？它就是把在循环中执行多次但结果却固定不变的代码移动到循环之外，这是最简单的优化。

第二个方法，减少过程调用。为什么要减少函数调用呢？因为每次调用一个函数，这个函数都可能会执行冗余的操作，比如数组边界检查，类型检查等。举例来说，如果我们需要调用一个函数来获得数组中的一个元素，我们可以首先获得数组首地址，然后通过索引访问每个元素，而非调用一个函数，然后检查数组是否越界，然后才返回元素值。

第三种方法，消除不必要的存储器引用。也就是说，在循环体中不使用指针来存放结果，而使用临时变量。为什么呢？目的就是减少读、写内存地址的时间。


随着对优化的深入，研究汇编代码以及理解现代处理器是如何执行计算则变得很重要。下面继续了解几种优化技术。
第四种方法，循环展开 。什么是循环展开？其思想是在一次迭代中对多个数组元素进行访问和合并操作。那这样做有什么原因或者说带来了什么好处呢？受限于资源约束，每个周期内只有两个功能单元可以执行数据相关的操作，然而一个周期内可能有涉及到数据的多个操作。那么，这些操作就不能在同一周期内完成，必须暂停等待，降低性能。通过循环展开可以在每次迭代中执行更多的数据操作，避免一个周期内产生多于两个操作，从而降低循环开销。

理想情况下的整数加法调度

实际资源受限的情况下，整数加法的调度

优化后，三次循环展开的整数加法的调度
 
第五种方法，迭代分割 。什么是迭代分割呢？我们知道处理器的几个功能单元是流水线化的，这说明它们在一个操作完成之前可以开始另一个新的操作。但是，我们的代码没有这种能力。当执行当前操作时，处理器会处理其它事务（一般不会暂停），直到本操作完成，才开始下一个操作。这严重影响了程序的性能。而通过迭代分割，则可以很好地缓解此问题，以提高性能。迭代分割其实就是把一个操作分割成两个或更多的部分，并在最后合并结果以提高性能。产生二路并行或n路并行，提高了程序的并行性。
 

理想情况下的整数乘法调度

实际资源受限的情况下，整数乘法的调度

优化后，二次展开、二次并行的整数乘法操作调度
比较直观的例子
for(;;){
  x0 = x0 + data[i];
  x1 = x1 + data[i+1];
}
result = x0 + x1;
运用上面的方法来优化第一部分多项式代码，思考！
int poly_eval(int *a, int degree, int x)
{
	int result = 0;
	int i;
	int xpwr = 1; /* Successive powers of x */
	for (i = 0; i <= degree; i++) {
		result += a[i]*xpwr;
		xpwr *= x;
	}
	return result;
}
对于第二部分的矩阵优化，我们用到的核心思想就是利用局部性原理，包括时间局部性和空间局部性，避免大幅度地跨越内存块。具体的示例在第三章“程序的机器级表示”的嵌套数组以及第六章“存储器层次结构”的局部性原理中讲述。
#define N 512
/* N x N matrix */
typedef int matrix_t[N][N];
/* Vector of length N */
typedef int vector_t[N];
void c_sum(matrix_t M, vector_t rowsum, vector_t colsum)
{
	int i,j;
	for (j = 0; j < N; j++) {
	colsum[j] = 0;
	for (i = 0; i < N; i++)
	colsum[j] += M[i][j];
	}
}
void rc_sum(matrix_t M, vector_t rowsum, vector_t colsum)
{
	int i,j;
	for (i = 0; i < N; i++) {
		rowsum[i] = colsum[i] = 0;
		for (j = 0; j < N; j++) {
		rowsum[i] += M[i][j];
		colsum[i] += M[j][i];
		}
	}
}