# MMA算法的推导及3D简支梁拓扑优化代码详解

5 篇文章 17 订阅

## MMA算法的推导及代码详解

### 问题描述

x = ( x 1 , x 2 , ⋅ ⋅ ⋅ x n ) T \mathbf{x} = \left( x_{1},x_{2}, \cdot \cdot \cdot x_{n} \right)^{T} 为n维欧氏空间 R n \mathbf{R}^{n} 内一点， f ( x ) , g i ( x ) ( i = 1 , 2 , ⋅ ⋅ ⋅ , m ) f\left( \mathbf{x} \right),g_{i}\left( \mathbf{x} \right)\left( i = 1,2, \cdot \cdot \cdot ,m \right) h j ( x ) ( j = m + 1 , ⋅ ⋅ ⋅ , p ) h_{j}\left( \mathbf{x} \right)\left( j = m + 1, \cdot \cdot \cdot ,p \right) 均为定的n元函数。最优化问题可表述为：在如下的约束条件下：
g i ( x ) ≤ 0      i = 1 , 2 ⋅ ⋅ ⋅ m g_{i}\left( \mathbf{x} \right) \leq 0~~~~i = 1,2 \cdot \cdot \cdot m

h j ( x ) = 0      j = m + 1 , ⋅ ⋅ ⋅ p h_{j}\left( \mathbf{x} \right) = 0~~~~j = m + 1, \cdot \cdot \cdot p

f ( x ) f(x) 的最小值（或最大值）。通常称 f ( x ) f(x) 为目标函数， g i ( x ) ≤ 0 g_{i}\left( \mathbf{x} \right) \leq 0 为不等式约束条件， h j ( x ) = 0 h_{j}\left( \mathbf{x} \right) = 0 为等式约束条件， x x 为设计变量。通用的表达形式为:

{ min ⁡      f ( x ) s . t .      g i ( x ) ≤ 0      i = 1 , 2 , ⋅ ⋅ ⋅ m   h j ( x ) = 0      j = m + 1 , ⋅ ⋅ ⋅ p x ∈ R n \left\{ \begin{matrix} {\min{~~~~f\left( \mathbf{x} \right)}} \\ {s.t.~~~~g_{i}\left( \mathbf{x} \right. )\leq 0~~~~i = 1,2, \cdot \cdot \cdot m~} \\ {h_{j}\left(\mathbf{x} \right.) = 0~~~~j = m + 1, \cdot \cdot \cdot p} \\ {x \in R^{n}} \\ \end{matrix} \right.

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\} x ∈ R x \in R 为上述问题的一个可行解集。

### 算法推导

P : P:
{ m i n       f 0 ( x ) f i ( x ) ≤ f i ^       f o r   i = 1 , … , m x j m i n ≤   x j ≤   x j m a x \begin{cases} min~~~~~f_0(x)\\ f_i(x)\le\widehat{f_i}~~~~~for~i=1,\dots,m\\ x_jmin\le~x_j\le~x_jmax \end{cases}

1. 选择合适的初始值
2. 计算迭代后得到的 x ( k ) x^{(k)} 以及每个约束函数的梯度
3. 产生原问题严格凸子问题 P ( k ) P^{(k)}
4. 利用迭代产生的结果求解子问题,并返回第二步

L j ( k ) < x j ( k ) < U j ( k ) L_j^{(k)}<x_j^{(k)}<U_j^{(k)}

f i ( k ) = r i ( k ) + ∑ j = 1 n ( p i j ( k ) U j ( k ) − x j + q i j ( k ) x j − L j ( k ) ) f_i^{(k)}=r_i^{(k)}+\sum_{j=1}^{n}(\frac{p_{ij}^{(k)}}{U_j^{(k)}-x_j}+\frac{q_{ij}^{(k)}}{x_j-L_j^{(k)}})
w h e r e     p i j ( k ) = { ( U j ( k ) − x j ( k ) ) 2 ∂ f i ∂ x j ,    i f ∂ f i x j > 0 0 ,                               i f ∂ f i x j ≤ 0 where~~~p_{ij}^{(k)}=\begin{cases} (U_{j}^{(k)}-x_j^{(k)})^2 \frac{\partial{f_i}}{\partial{x_j}},~~if\frac{\partial{f_i}}{x_j}>0\\ 0,~~~~~~~~~~~~~~~~~~~~~~~~~~~~~if\frac{\partial{f_i}}{x_j}\le0 \end{cases}
q i j ( k ) = { 0 ,                                i f ∂ f i x j ≥ 0 − ( x j ( k ) − L j ( k ) ) 2 ∂ f i ∂ x j ,    i f ∂ f i x j < 0 q_{ij}^{(k)}=\begin{cases} 0,~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~if\frac{\partial{f_i}}{x_j}\geq0\\ -(x_{j}^{(k)}-L_j^{(k)})^2 \frac{\partial{f_i}}{\partial{x_j}},~~if\frac{\partial{f_i}}{x_j}<0\\ \end{cases}

r i ( k ) = f i ( x ( k ) ) − ∑ j = 1 n ( p i j ( k ) U j ( k ) − x j ( k ) + q i j ( k ) x j ( k ) − L j ( k ) ) r_i^{(k)}=f_i(x^{(k)})-\sum_{j=1}^{n}(\frac{p_{ij}^{(k)}}{U_j^{(k)}-x_j^{(k)}}+\frac{q_{ij}^{(k)}}{x_j^{(k)}-L_j^{(k)}})

f i ( k ) = f i ( x ( k ) ) + ∂ f i ∂ x ⋅ ( x − x k ) ⋅ ( U − x k ) ( U − x ) f_i^{(k)}=f_i(x^{(k)})+\frac{\partial{f_i}}{\partial{x}}\cdot{(x-x^{k})}\cdot\frac{(U-x^k)}{(U-x)}
= f i ( x ( k ) ) + ∂ f i ∂ x ⋅ ( x − x k ) =f_i(x^{(k)})+\frac{\partial{f_i}}{\partial{x}}\cdot{(x-x^{k})}

L → − ∞ L\to{-\infty} 可以得到相似的结论。

∂ 2 f i ( k ) ∂ x j ( 2 ) = 2 p i j ( k ) ( U j ( k ) − x j ) 3 + 2 q i j ( k ) ( x j − L j ( k ) ) 3 \frac{\partial^2{f_i^{(k)}}}{\partial{x_j^{(2)}}}=\frac{2p_{ij^{(k)}}}{(U_j^{(k)}-x_j)^3}+\frac{2q_{ij}^{(k)}}{(x_j-L_j^{(k)})^{3}}

∂ 2 f i ( k ) ∂ x j ∂ x i = 0 , i f   i ≠ j \frac{\partial^2f_i^{(k)}}{\partial{x_j}\partial{x_i}}=0,if~i\neq{j}

L j ( k ) = x j ( k ) − ( x j m a x − x j m i n )    &    U j ( k ) = x j ( k ) + ( x j m a x − x j m i n ) L^{(k)}_j=x_j^{(k)}-(x_jmax-x_jmin)~~\&~~U^{(k)}_j=x_j^{(k)}+(x_jmax-x_jmin)

L j ( k ) = x j ( k ) − s ( x j ( k − 1 ) − L j ( k − 1 ) )    &    U j ( k ) = x j ( k ) + s ( x j ( k − 1 ) − U j ( k − 1 ) ) L^{(k)}_j=x_j^{(k)}-s(x_j^{(k-1)}-L_j^{(k-1)})~~\&~~U^{(k)}_j=x_j^{(k)}+s(x_j^{(k-1)}-U_j^{(k-1)})

m i n                   ∑ j = 1 n ( p 0 j U j − x j − q 0 j x j − L j ) + r 0 min~~~~~~~~~~~~~~~~~\sum_{j=1}^{n}(\frac{p_0j}{U_j-x_j}-\frac{q_0j}{x_j-L_j})+r_0
s . t . ∑ j = 1 n ( p i j U j − x j − q i j x j − L j ) ≤ b i s.t.\sum_{j=1}^{n}(\frac{p_{ij}}{U_j-x_j}-\frac{q_{ij}}{x_j-L_j})\le{b_i}
α j ≤ x j ≤ β j \alpha_{j}\leq{x_j}\leq{\beta_j}

l ( x , y ) = f 0 ( k ) ( x ) + ∑ j = 1 n y i f i ( k ) ( x ) l(x,y)=f_0^{(k)}(x)+\sum_{j=1}^{n}y_if_i^{(k)}(x)

l ( x , y ) = p 0 j + y T p j U j − x j − q 0 j + y T q j x j − L j l(x,y)=\frac{p_{0j}+y^Tp_j}{U_j-x_j}-\frac{q_{0j}+y^Tq_j}{x_j-L_j}

W ( y ) = m i n x { L ( x , y ) ; α j ≤ x j ≤ β j , f o r   a l l   j } W(y)=\underset{x}{min}\{{L(x,y);\alpha_{j}\leq{x_j}\leq{\beta_j}},for~all~j\}

P ( K ) P^{(K)} 可以转化为对偶问题:
m a x      W ( y ) = r 0 − y T b + ∑ j = 1 n ( p 0 j + y T p j U j − x j − q 0 j + y T q j x j − L j ) max~~~~W(y)=r_0-y^Tb+\sum_{j=1}^{n}(\frac{p_{0j}+y^Tp_j}{U_j-x_j}-\frac{q_{0j}+y^Tq_j}{x_j-L_j})
s . t .         y ≥ 0 s.t.~~~~~~~y\geq{0}

m i n i m i z e            ∑ j = 1 n ( p 0 j U j − x j − q 0 j x j − L j ) + ∑ i = 1 m ( d i z i + d i z i 2 ) + r 0 minimize~~~~~~~~~~\sum_{j=1}^{n}(\frac{p_{0j}}{U_j-x_j}-\frac{q_{0j}}{x_j-L_j})+\sum_{i=1}^{m}(d_iz_i+d_iz_i^2)+r_0
s . t .                      ∑ j = 1 n ( p i j U j − x j − q i j x j − L j ) − z i ≤ b i ,   f o r   i = 1 , . . . , m s.t. ~~~~~~~~~~~~~~~~~~~~\sum_{j=1}^{n}(\frac{p_{ij}}{U_j-x_j}-\frac{q_{ij}}{x_j-L_j})-z_i\leq{b_i},~for~i=1,...,m
a n d                      α j ≤ x j ≤ β j ,   f o r   j = 1 , . . . , n   , z i ≥ 0 ,   f o r   i = 1 , . . . , m and~~~~~~~~~~~~~~~~~~~~\alpha_{j}\leq{x_j}\leq{\beta_j},~for~j=1,...,n~,z_i\ge{0},~for~i=1,...,m

### 代码及注释

top3d

%本文给定top3d(30,20,10,0.4,3,1.5)
function top3d(nelx,nely,nelz,volfrac,penal,rmin)
% 定义循环参数
maxloop = 200;    % 最大迭代次数
tolx = 0.01;      % 终止条件（残差）
displayflag = 0;  % 显示结构表示
% 材料的属性
E0 = 1;           % 固体区域的杨氏模量
Emin = 1e-9;      % 非固体区域的杨氏模量，尽可能小但是为了避免奇异性一般不取0
nu = 0.3;         % 泊松比
[il,jl,kl] = meshgrid(nelx, 0, 0:nelz);                 % Coordinates
loadnid = kl*(nelx+1)*(nely+1)+il*(nely+1)+(nely+1-jl); % Node IDs
% 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);
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;
% 开始迭代
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 %#ok<UNRCH>
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


mma

%-------------------------------------------------------
% 	 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个α<x_j对应的拉格朗日乘子
%  eta   = n个x_j<β对应的拉格朗日乘子
%   mu   = m个y_i<0对应的拉格朗日乘子。
%  zet   = z>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);


subsolve

%  子问题求解程序 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;


### 参考文献

[1]SVANBERG K. MMA and GCMMA – two methods for nonlinear optimization[J]. : 15.

[2]STEPHEN BOYD, LIEVEN VANDENBERGHE. 凸优化[M]. 王书宁, 译, 许鋆, 译, 黄晓霖, 译. 清华大学出版社.

[3]LIU K, TOVAR A. An efficient 3D topology optimization code written in Matlab[J]. Structural and Multidisciplinary Optimization, 2014, 50(6): 1175–1196.

• 14
点赞
• 75
收藏
觉得还不错? 一键收藏
• 15
评论
03-06
01-08
05-03 1万+
12-18 2602
06-01 3003
07-28
09-30

### “相关推荐”对你有帮助么？

• 非常没帮助
• 没帮助
• 一般
• 有帮助
• 非常有帮助

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。