diff --git a/MMA_test/ISOdisplay.fig b/MMA_test/ISOdisplay.fig new file mode 100644 index 0000000..005e536 Binary files /dev/null and b/MMA_test/ISOdisplay.fig differ diff --git a/MMA_test/mmasub.m b/MMA_test/mmasub.m new file mode 100644 index 0000000..a98c215 --- /dev/null +++ b/MMA_test/mmasub.m @@ -0,0 +1,133 @@ +%------------------------------------------------------- +% MMA求解程序 MMA.m +% 本代码将原本的非线性问题线性近似，非凸问题凸近似，得到一系列 +% 线性的、凸的子问题，再通过subsolv.m迭代求解这些线性的、凸的子问题最终找到最优值点 +function [xmma,ymma,zmma,lam,xsi,eta,mu,zet,s,low,upp] = ... +mmasub(m,n,iter,xval,xmin,xmax,xold1,xold2, ... +f0val,df0dx,fval,dfdx,low,upp,a0,a,c,d) +% This function mmasub performs one MMA-iteration, aimed at +% solving the nonlinear programming problem: +% 该函数进行MMA迭代，已解决符合以下形式的非线性（非凸）优化问题 +% 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 +% 输入参数有： +% m = 约束条件的数量 +% n = 变量的个数 +% iter = 迭代步数 +% xval = 包含所有xj的向量 +% xmin = 包含所有xj取值下限的向量 +% xmax = 包含所有xj取值上限的向量 +% xold1 = 上次迭代产生的x_val +% xold2 = 上上次迭代产生的的xval +% f0val = 目标函数的现有值 +% df0dx = 包含目标函数对所有xj偏导数的列向量 +% fval = 时约束函数fi值的列向量 +% dfdx = x_j.fi对xj偏导数的矩阵 +% low = 前一次迭代的渐近线下界 +% upp = 前一次迭代的渐进线上界 +% 引入下面四个量防止在初始值选择不恰当的形况下出现子问题不适定（没有结）的情况 +% a0 = a_0*z中的常系数a_0 +% a = 用来存储a_i*z中系数 a_i的列向量 +% c = 用来存储c_i*y_i中系数c_i的列向量 +% d =用来存储0.5*d_i*(y_i)^2中系数d_i的列向量，d_i应该是一个足够大的数， +% 可以证明当d_i足够大时y_i趋近于0 +%输出的量: +% xmma = 子问题求解出的xj的列向量 +% ymma = 子问题求解出的yi的列向量 +% zmma = z的标量 +% lam = m个约束对应的拉格朗日乘子 +% xsi = n个α0对应的拉格朗日乘子 +% s = m个约束条件对应的松弛因子 +% low = 本次求解的子问题对应的渐近线的下界 +% upp = 本次求解的子问题对应的渐近线上界 +% +%epsimin = sqrt(m+n)*10^(-9);计算机精度的最小值 +epsimin = 10^(-7); +raa0 = 0.00001; +move = 1.0; +albefa = 0.1; +asyinit = 0.5; +asyincr = 1.2; +asydecr = 0.7; +eeen = ones(n,1); +eeem = ones(m,1); +zeron = zeros(n,1); + +% Calculation of the asymptotes low and upp :计算渐近线的上下界 +if iter < 2.5 + low = xval - asyinit*(xmax-xmin); + upp = xval + asyinit*(xmax-xmin); +else + zzz = (xval-xold1).*(xold1-xold2); + factor = eeen; + factor(find(zzz > 0)) = asyincr; + factor(find(zzz < 0)) = asydecr; + low = xval - factor.*(xold1 - low); + upp = xval + factor.*(upp - xold1); + lowmin = xval - 10*(xmax-xmin); + lowmax = xval - 0.01*(xmax-xmin); + uppmin = xval + 0.01*(xmax-xmin); + uppmax = xval + 10*(xmax-xmin); + low = max(low,lowmin); + low = min(low,lowmax); + upp = min(upp,uppmax); + upp = max(upp,uppmin); +end + +% 计算设计变量的精确界限α和β + +zzz1 = low + albefa*(xval-low); +zzz2 = xval - move*(xmax-xmin); +zzz = max(zzz1,zzz2); +alfa = max(zzz,xmin); +zzz1 = upp - albefa*(upp-xval); +zzz2 = xval + move*(xmax-xmin); +zzz = min(zzz1,zzz2); +beta = min(zzz,xmax); + +% 计算 p0, q0, P, Q 和 b. + +xmami = xmax-xmin; +xmamieps = 0.00001*eeen; +xmami = max(xmami,xmamieps); +xmamiinv = eeen./xmami; +ux1 = upp-xval; +ux2 = ux1.*ux1; +xl1 = xval-low; +xl2 = xl1.*xl1; +uxinv = eeen./ux1; +xlinv = eeen./xl1; +% +p0 = zeron; +q0 = zeron; +p0 = max(df0dx,0); +q0 = max(-df0dx,0); +%p0(find(df0dx > 0)) = df0dx(find(df0dx > 0)); +%q0(find(df0dx < 0)) = -df0dx(find(df0dx < 0)); +pq0 = 0.001*(p0 + q0) + raa0*xmamiinv; +p0 = p0 + pq0; +q0 = q0 + pq0; +p0 = p0.*ux2; +q0 = q0.*xl2; +% +P = sparse(m,n); +Q = sparse(m,n); +P = max(dfdx,0); +Q = max(-dfdx,0); +%P(find(dfdx > 0)) = dfdx(find(dfdx > 0)); +%Q(find(dfdx < 0)) = -dfdx(find(dfdx < 0)); +PQ = 0.001*(P + Q) + raa0*eeem*xmamiinv'; +P = P + PQ; +Q = Q + PQ; +P = P * spdiags(ux2,0,n,n); +Q = Q * spdiags(xl2,0,n,n); +b = P*uxinv + Q*xlinv - fval ; +% +%%% 使用反演原对偶牛顿方法求解子问题 +[xmma,ymma,zmma,lam,xsi,eta,mu,zet,s] = ... +subsolv(m,n,epsimin,low,upp,alfa,beta,p0,q0,P,Q,a0,a,b,c,d); diff --git a/MMA_test/subsolv.m b/MMA_test/subsolv.m new file mode 100644 index 0000000..c296e72 --- /dev/null +++ b/MMA_test/subsolv.m @@ -0,0 +1,211 @@ +% 子问题求解程序 subsolv.m +% +function [xmma,ymma,zmma,lamma,xsimma,etamma,mumma,zetmma,smma] = ... +subsolv(m,n,epsimin,low,upp,alfa,beta,p0,q0,P,Q,a0,a,b,c,d); +% +% 使用该函数求解MMA生成的子问题可以写成： +% +% minimize SUM[ p0j/(uppj-xj) + q0j/(xj-lowj) ] + a0*z + +% + SUM[ ci*yi + 0.5*di*(yi)^2 ], +% +% subject to SUM[ pij/(uppj-xj) + qij/(xj-lowj) ] - ai*z - yi <= bi, +% alfaj <= xj <= betaj, yi >= 0, z >= 0. +% 然后使用拉格朗日的方法求解子问题的对偶问题 +% Input: m, n, low, upp, alfa, beta, p0, q0, P, Q, a0, a, b, c, d. +% Output: xmma,ymma,zmma, slack variables and Lagrange multiplers. +% +een = ones(n,1); +eem = ones(m,1); +epsi = 1; +epsvecn = epsi*een; +epsvecm = epsi*eem; +x = 0.5*(alfa+beta); +y = eem; +z = 1; +lam = eem; +xsi = een./(x-alfa); +xsi = max(xsi,een); +eta = een./(beta-x); +eta = max(eta,een); +mu = max(eem,0.5*c); +zet = 1; +s = eem; +itera = 0; + +while epsi > epsimin + epsvecn = epsi*een; + epsvecm = epsi*eem; + ux1 = upp-x; + xl1 = x-low; + ux2 = ux1.*ux1; + xl2 = xl1.*xl1; + uxinv1 = een./ux1; + xlinv1 = een./xl1; + + plam = p0 + P'*lam ; + qlam = q0 + Q'*lam ; + gvec = P*uxinv1 + Q*xlinv1; + dpsidx = plam./ux2 - qlam./xl2 ; + + rex = dpsidx - xsi + eta; + rey = c + d.*y - mu - lam; + rez = a0 - zet - a'*lam; + relam = gvec - a*z - y + s - b; + rexsi = xsi.*(x-alfa) - epsvecn; + reeta = eta.*(beta-x) - epsvecn; + remu = mu.*y - epsvecm; + rezet = zet*z - epsi; + res = lam.*s - epsvecm; + + residu1 = [rex' rey' rez]'; + residu2 = [relam' rexsi' reeta' remu' rezet res']'; + residu = [residu1' residu2']'; + residunorm = sqrt(residu'*residu); + residumax = max(abs(residu)); + + ittt = 0; + while residumax > 0.9*epsi & ittt < 100 + ittt=ittt + 1; + itera=itera + 1; + + ux1 = upp-x; + xl1 = x-low; + ux2 = ux1.*ux1; + xl2 = xl1.*xl1; + ux3 = ux1.*ux2; + xl3 = xl1.*xl2; + uxinv1 = een./ux1; + xlinv1 = een./xl1; + uxinv2 = een./ux2; + xlinv2 = een./xl2; + plam = p0 + P'*lam ; + qlam = q0 + Q'*lam ; + gvec = P*uxinv1 + Q*xlinv1; + GG = P*spdiags(uxinv2,0,n,n) - Q*spdiags(xlinv2,0,n,n); + dpsidx = plam./ux2 - qlam./xl2 ; + delx = dpsidx - epsvecn./(x-alfa) + epsvecn./(beta-x); + dely = c + d.*y - lam - epsvecm./y; + delz = a0 - a'*lam - epsi/z; + dellam = gvec - a*z - y - b + epsvecm./lam; + diagx = plam./ux3 + qlam./xl3; + diagx = 2*diagx + xsi./(x-alfa) + eta./(beta-x); + diagxinv = een./diagx; + diagy = d + mu./y; + diagyinv = eem./diagy; + diaglam = s./lam; + diaglamyi = diaglam+diagyinv; + + if m < n + blam = dellam + dely./diagy - GG*(delx./diagx); + bb = [blam' delz]'; + Alam = spdiags(diaglamyi,0,m,m) + GG*spdiags(diagxinv,0,n,n)*GG'; + AA = [Alam a + a' -zet/z ]; + solut = AA\bb; + dlam = solut(1:m); + dz = solut(m+1); + dx = -delx./diagx - (GG'*dlam)./diagx; + else + diaglamyiinv = eem./diaglamyi; + dellamyi = dellam + dely./diagy; + Axx = spdiags(diagx,0,n,n) + GG'*spdiags(diaglamyiinv,0,m,m)*GG; + azz = zet/z + a'*(a./diaglamyi); + axz = -GG'*(a./diaglamyi); + bx = delx + GG'*(dellamyi./diaglamyi); + bz = delz - a'*(dellamyi./diaglamyi); + AA = [Axx axz + axz' azz ]; + bb = [-bx' -bz]'; + solut = AA\bb; + dx = solut(1:n); + dz = solut(n+1); + dlam = (GG*dx)./diaglamyi - dz*(a./diaglamyi) + dellamyi./diaglamyi; + end + + dy = -dely./diagy + dlam./diagy; + dxsi = -xsi + epsvecn./(x-alfa) - (xsi.*dx)./(x-alfa); + deta = -eta + epsvecn./(beta-x) + (eta.*dx)./(beta-x); + dmu = -mu + epsvecm./y - (mu.*dy)./y; + dzet = -zet + epsi/z - zet*dz/z; + ds = -s + epsvecm./lam - (s.*dlam)./lam; + xx = [ y' z lam' xsi' eta' mu' zet s']'; + dxx = [dy' dz dlam' dxsi' deta' dmu' dzet ds']'; + + stepxx = -1.01*dxx./xx; + stmxx = max(stepxx); + stepalfa = -1.01*dx./(x-alfa); + stmalfa = max(stepalfa); + stepbeta = 1.01*dx./(beta-x); + stmbeta = max(stepbeta); + stmalbe = max(stmalfa,stmbeta); + stmalbexx = max(stmalbe,stmxx); + stminv = max(stmalbexx,1); + steg = 1/stminv; + + xold = x; + yold = y; + zold = z; + lamold = lam; + xsiold = xsi; + etaold = eta; + muold = mu; + zetold = zet; + sold = s; + + itto = 0; + resinew = 2*residunorm; + while resinew > residunorm & itto < 50 + itto = itto+1; + + x = xold + steg*dx; + y = yold + steg*dy; + z = zold + steg*dz; + lam = lamold + steg*dlam; + xsi = xsiold + steg*dxsi; + eta = etaold + steg*deta; + mu = muold + steg*dmu; + zet = zetold + steg*dzet; + s = sold + steg*ds; + ux1 = upp-x; + xl1 = x-low; + ux2 = ux1.*ux1; + xl2 = xl1.*xl1; + uxinv1 = een./ux1; + xlinv1 = een./xl1; + plam = p0 + P'*lam ; + qlam = q0 + Q'*lam ; + gvec = P*uxinv1 + Q*xlinv1; + dpsidx = plam./ux2 - qlam./xl2 ; + + rex = dpsidx - xsi + eta; + rey = c + d.*y - mu - lam; + rez = a0 - zet - a'*lam; + relam = gvec - a*z - y + s - b; + rexsi = xsi.*(x-alfa) - epsvecn; + reeta = eta.*(beta-x) - epsvecn; + remu = mu.*y - epsvecm; + rezet = zet*z - epsi; + res = lam.*s - epsvecm; + + residu1 = [rex' rey' rez]'; + residu2 = [relam' rexsi' reeta' remu' rezet res']'; + residu = [residu1' residu2']'; + resinew = sqrt(residu'*residu); + steg = steg/2; + end + residunorm=resinew; + residumax = max(abs(residu)); + steg = 2*steg; + end +epsi = 0.1*epsi; +end + +xmma = x; +ymma = y; +zmma = z; +lamma = lam; +xsimma = xsi; +etamma = eta; +mumma = mu; +zetmma = zet; +smma = s; diff --git a/MMA_test/top3d.m b/MMA_test/top3d.m new file mode 100644 index 0000000..04afd35 --- /dev/null +++ b/MMA_test/top3d.m @@ -0,0 +1,190 @@ +%本文给定top3d(30,20,10,0.4,3,1.5) +function top3d(nelx,nely,nelz,volfrac,penal,rmin) +% 定义循环参数 +maxloop = 200; % 最大迭代次数 +tolx = 0.01; % 终止条件（残差） +displayflag = 1; % 显示结构表示 +% 材料的属性 +E0 = 1; % 固体区域的杨氏模量 +Emin = 1e-9; % 非固体区域的杨氏模量，尽可能小但是为了避免奇异性一般不取0 +nu = 0.3; % 泊松比 +% USER-DEFINED LOAD DOFs +[il,jl,kl] = meshgrid(nelx, 0, 0:nelz); % Coordinates +loadnid = kl*(nelx+1)*(nely+1)+il*(nely+1)+(nely+1-jl); % Node IDs +loaddof = 3*loadnid(:) - 1; % DOFs +% USER-DEFINED SUPPORT FIXED DOFs +[iif,jf,kf] = meshgrid(0,0:nely,0:nelz); % Coordinates +fixednid = kf*(nelx+1)*(nely+1)+iif*(nely+1)+(nely+1-jf); % Node IDs +fixeddof = [3*fixednid(:); 3*fixednid(:)-1; 3*fixednid(:)-2]; % DOFs +% 有限元分析程序 +nele = nelx*nely*nelz; +ndof = 3*(nelx+1)*(nely+1)*(nelz+1); +F = sparse(loaddof,1,-1,ndof,1); +U = zeros(ndof,1); +freedofs = setdiff(1:ndof,fixeddof); +KE = lk_H8(nu); +nodegrd = reshape(1:(nely+1)*(nelx+1),nely+1,nelx+1); +nodeids = reshape(nodegrd(1:end-1,1:end-1),nely*nelx,1); +nodeidz = 0:(nely+1)*(nelx+1):(nelz-1)*(nely+1)*(nelx+1); +nodeids = repmat(nodeids,size(nodeidz))+repmat(nodeidz,size(nodeids)); +edofVec = 3*nodeids(:)+1; +edofMat = repmat(edofVec,1,24)+ ... + repmat([0 1 2 3*nely + [3 4 5 0 1 2] -3 -2 -1 ... + 3*(nely+1)*(nelx+1)+[0 1 2 3*nely + [3 4 5 0 1 2] -3 -2 -1]],nele,1); +iK = reshape(kron(edofMat,ones(24,1))',24*24*nele,1); +jK = reshape(kron(edofMat,ones(1,24))',24*24*nele,1); +% 过滤器 +iH = ones(nele*(2*(ceil(rmin)-1)+1)^2,1); +jH = ones(size(iH)); +sH = zeros(size(iH)); +k = 0; +for k1 = 1:nelz + for i1 = 1:nelx + for j1 = 1:nely + e1 = (k1-1)*nelx*nely + (i1-1)*nely+j1; + for k2 = max(k1-(ceil(rmin)-1),1):min(k1+(ceil(rmin)-1),nelz) + 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 = (k2-1)*nelx*nely + (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+(k1-k2)^2)); + end + end + end + end + end +end +H = sparse(iH,jH,sH); +Hs = sum(H,2); +% 迭代初始化 +x = repmat(volfrac,[nely,nelx,nelz]); +xPhys = x; +loop = 0; +change = 1; + +m=1; +n=nele; +xmin=zeros (n, 1) ; +xmax=ones (n, 1) ; +xold1=x(:); +xold2=x(:); +low=ones(n, 1) ; +upp=ones(n, 1) ; +a0=1; +a=zeros (m, 1) ; +c_MMA=10000*ones (m, 1); +d=zeros (m, 1) ; +% 开始迭代 +while change > tolx && loop < maxloop + loop = loop+1; + % 有限元分析 + sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(E0-Emin)),24*24*nele,1); + K = sparse(iK,jK,sK); K = (K+K')/2; + U(freedofs,:) = K(freedofs,freedofs)\F(freedofs,:); + % 定义目标函数与灵敏度分析 + ce = reshape(sum((U(edofMat)*KE).*U(edofMat),2),[nely,nelx,nelz]); + c = sum(sum(sum((Emin+xPhys.^penal*(E0-Emin)).*ce))); + dc = -penal*(E0-Emin)*xPhys.^(penal-1).*ce; + dv = ones(nely,nelx,nelz); + % 对灵敏度进行过滤与投影，需要注意的是，并不能数学上证明对灵敏度进行过滤的 + % 稳定性，如果出现不适定的情况可以尝试更改过滤为密度过滤 + dc(:) = H*(dc(:)./Hs); + dv(:) = H*(dv(:)./Hs); + % 移动渐近线法求解 + xval = x(:); + f0val = c; + df0dx = dc(:); + fval = sum(xPhys(:))/(volfrac*nele) - 1; + dfdx = dv(:)' / (volfrac*nele); + [xmma, ~, ~, ~, ~, ~, ~, ~, ~, low,upp] = ... + mmasub(m, n, loop, xval, xmin, xmax, xold1, xold2, ... + f0val,df0dx,fval,dfdx,low,upp,a0,a,c_MMA,d); + % 更新迭代变量 + xnew = reshape(xmma,nely,nelx,nelz); + xPhys(:) = (H*xnew(:))./Hs; + xold2 = xold1(:); + xold1 = x(:); + change = max(abs(xnew(:)-x(:))); + x = xnew; + % 输出结果 + fprintf(' It.:%5i Obj.:%11.4f Vol.:%7.3f ch.:%7.3f\n',loop,c,mean(xPhys(:)),change); + % 输出密度 + if displayflag, clf; display_3D(xPhys); end +end +clf; display_3D(xPhys); +end + + +% === 生成单元刚度矩阵 === +function [KE] = lk_H8(nu) +A = [32 6 -8 6 -6 4 3 -6 -10 3 -3 -3 -4 -8; + -48 0 0 -24 24 0 0 0 12 -12 0 12 12 12]; +k = 1/144*A'*[1; nu]; + +K1 = [k(1) k(2) k(2) k(3) k(5) k(5); + k(2) k(1) k(2) k(4) k(6) k(7); + k(2) k(2) k(1) k(4) k(7) k(6); + k(3) k(4) k(4) k(1) k(8) k(8); + k(5) k(6) k(7) k(8) k(1) k(2); + k(5) k(7) k(6) k(8) k(2) k(1)]; +K2 = [k(9) k(8) k(12) k(6) k(4) k(7); + k(8) k(9) k(12) k(5) k(3) k(5); + k(10) k(10) k(13) k(7) k(4) k(6); + k(6) k(5) k(11) k(9) k(2) k(10); + k(4) k(3) k(5) k(2) k(9) k(12) + k(11) k(4) k(6) k(12) k(10) k(13)]; +K3 = [k(6) k(7) k(4) k(9) k(12) k(8); + k(7) k(6) k(4) k(10) k(13) k(10); + k(5) k(5) k(3) k(8) k(12) k(9); + k(9) k(10) k(2) k(6) k(11) k(5); + k(12) k(13) k(10) k(11) k(6) k(4); + k(2) k(12) k(9) k(4) k(5) k(3)]; +K4 = [k(14) k(11) k(11) k(13) k(10) k(10); + k(11) k(14) k(11) k(12) k(9) k(8); + k(11) k(11) k(14) k(12) k(8) k(9); + k(13) k(12) k(12) k(14) k(7) k(7); + k(10) k(9) k(8) k(7) k(14) k(11); + k(10) k(8) k(9) k(7) k(11) k(14)]; +K5 = [k(1) k(2) k(8) k(3) k(5) k(4); + k(2) k(1) k(8) k(4) k(6) k(11); + k(8) k(8) k(1) k(5) k(11) k(6); + k(3) k(4) k(5) k(1) k(8) k(2); + k(5) k(6) k(11) k(8) k(1) k(8); + k(4) k(11) k(6) k(2) k(8) k(1)]; +K6 = [k(14) k(11) k(7) k(13) k(10) k(12); + k(11) k(14) k(7) k(12) k(9) k(2); + k(7) k(7) k(14) k(10) k(2) k(9); + k(13) k(12) k(10) k(14) k(7) k(11); + k(10) k(9) k(2) k(7) k(14) k(7); + k(12) k(2) k(9) k(11) k(7) k(14)]; +KE = 1/((nu+1)*(1-2*nu))*... + [ K1 K2 K3 K4; + K2' K5 K6 K3'; + K3' K6 K5' K2'; + K4 K3 K2 K1']; +end +% === 展示3D效果图 (ISO-VIEW) === +function display_3D(rho) +[nely,nelx,nelz] = size(rho); +hx = 1; hy = 1; hz = 1; % 定义效果图的单元大小 +face = [1 2 3 4; 2 6 7 3; 4 3 7 8; 1 5 8 4; 1 2 6 5; 5 6 7 8]; +set(gcf,'Name','ISO display','NumberTitle','off'); +for k = 1:nelz + z = (k-1)*hz; + for i = 1:nelx + x = (i-1)*hx; + for j = 1:nely + y = nely*hy - (j-1)*hy; + if (rho(j,i,k) > 0.5) % 定义展示出的密度阈值 + vert = [x y z; x y-hx z; x+hx y-hx z; x+hx y z; x y z+hx;x y-hx z+hx; x+hx y-hx z+hx;x+hx y z+hx]; + vert(:,[2 3]) = vert(:,[3 2]); vert(:,2,:) = -vert(:,2,:); + patch('Faces',face,'Vertices',vert,'FaceColor',[0.2+0.8*(1-rho(j,i,k)),0.2+0.8*(1-rho(j,i,k)),0.2+0.8*(1-rho(j,i,k))]); + hold on; + end + end + end +end +axis equal; axis tight; axis off; box on; view([30,30]); pause(1e-6); +end