拓扑优化 mma法

%=========================================================================%
% 2019.03.15 BIT 宇航学院 DONGJZ
% 参考：
% [1]https://www.researchgate.net/profile/Ole_Sigmund
% [3]博客园
%=========================================================================%
clc;clear all;format long g;warning off all;
global volfrac nelx nely h1 h2 nyi nyj nxi nxj NFORCE NFIXED
global YOUNG POSSION THICK
global FORCE FIXED
%=========================================================================%
% 变量名称
%=========================================================================%
% volfrac-质量约束，nelx*nely=单元个数，h1*h2=网格面积
% NFORCE -受力节点总数，NFIXED-固定节点总数
% YOUUNG-模量，POISS-泊松比，THICK-厚度
% nyi nyj nxi nxj 空洞单元区域，i为下届，j为上届 %此段有问题 初置零；空洞区域
% 有限元程序没有搞出来，不均匀场太大，刚度阵奇异了？？？怎么办？
% FORCE-节点力数组（n*3）={受力节点，x方向，y方向}
% FIXED-约束信息数组（n*3）={约束点，x方向，y方向}1为约束，否则0
% FP1 数据文件指针
% 拓扑优化很简单，大部分都已经软件化了，这一块的创新在于算法，可以说第一篇关于
%拓扑优化设计的思想，就好比刘正猷的mass in mass模型，是开创篇，Sigmund确实厉害，
%能够整出个MMA不一般啊！！！这些个东西都很简单，基本的数学知识，现在想来，初高中
%的数学教育是有点用的。解决问题，好像真的是最简单的就是不错的。
%在oc求解器的基准上，做了个mma嫁接。mama对非线性函数有很大的帮助，对于线性约束，
%先转化为非线性约束，这个好像是复杂化了，不过只是一个探讨。
%=========================================================================%
% 输入相关数据
%=========================================================================%
FP1=fopen('input.txt','rt');
volfrac=fscanf(FP1,'%f',1);
nely=fscanf(FP1,'%f',1);
nelx=fscanf(FP1,'%f',1);
h1=fscanf(FP1,'%f',1);
h2=fscanf(FP1,'%f',1);
nyi=fscanf(FP1,'%f',1);nyj=fscanf(FP1,'%f',1);
nxi=fscanf(FP1,'%f',1);nxj=fscanf(FP1,'%f',1);
NFORCE=fscanf(FP1,'%d',1);NFIXED=fscanf(FP1,'%d',1);
YOUNG=fscanf(FP1,'%f',1);
POSSION=fscanf(FP1,'%f',1);
THICK=fscanf(FP1,'%f',1);
FORCE=fscanf(FP1,'%f',[3,NFORCE])';%节点力数组
FIXED=fscanf(FP1,'%f',[3,NFIXED])';%约束信息
%==========================================================================%
% 初始化
%==========================================================================%
penal=5;
rmin=h1*1.2;
%==========================================================================%
% Minimize f_ 0(x) + a_ 0*z + sum( c_i*y_i + 0.5*d_i*(y_i)^2 )
% subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m
% xmin_j <= x_j <= xmax_j, j = 1,...,n
% z >= 0, y_i >= 0, i = 1,...,m
%==========================================================================%
nx=nely*nelx;nm=1;
xcur(1:nely,1:nelx)=0.3;
xold1=volfrac*ones(nx,1);
xold2=volfrac*ones(nx,1);
xmin=0.001*ones(nx,1);xmax=ones(nx,1);
lowcur=zeros(nx,1);uppcur=zeros(nx,1);
a0=0;a=0;c=0;d=0;
%==========================================================================%
% 开始迭代
%=====&