RBTO-PMA-DLM_RBTO_拓扑优化_可靠度优化_

%------------------------------------------------------- % This is the file mmasub.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); % % Version September 2007 (and a small change August 2008) % % Krister Svanberg <krille@math.kth.se> % Department of Mathematics, SE-10044 Stockholm, Sweden. % % This function mmasub performs one MMA-iteration, aimed at % solving the nonlinear programming problem: % % 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 %*** INPUT: % % m = The number of general constraints. % n = The number of variables x_j. % iter = Current iteration number ( =1 the first time mmasub is called). % xval = Column vector with the current values of the variables x_j. % xmin = Column vector with the lower bounds for the variables x_j. % xmax = Column vector with the upper bounds for the variables x_j. % xold1 = xval, one iteration ago (provided that iter>1). % xold2 = xval, two iterations ago (provided that iter>2). % f0val = The value of the objective function f_0 at xval. % df0dx = Column vector with the derivatives of the objective function % f_0 with respect to the variables x_j, calculated at xval. % fval = Column vector with the values of the constraint functions f_i, % calculated at xval. % dfdx = (m x n)-matrix with the derivatives of the constraint functions % f_i with respect to the variables x_j, calculated at xval. % dfdx(i,j) = the derivative of f_i with respect to x_j. % low = Column vector with the lower asymptotes from the previous % iteration (provided that iter>1). % upp = Column vector with the upper asymptotes from the previous % iteration (provided that iter>1). % a0 = The constants a_0 in the term a_0*z. % a = Column vector with the constants a_i in the terms a_i*z. % c = Column vector with the constants c_i in the terms c_i*y_i. % d = Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2. % %*** OUTPUT: % % xmma = Column vector with the optimal values of the variables x_j % in the current MMA subproblem. % ymma = Column vector with the optimal values of the variables y_i % in the current MMA subproblem. % zmma = Scalar with the optimal value of the variable z % in the current MMA subproblem. % lam = Lagrange multipliers for the m general MMA constraints. % xsi = Lagrange multipliers for the n constraints alfa_j - x_j <= 0. % eta = Lagrange multipliers for the n constraints x_j - beta_j <= 0. % mu = Lagrange multipliers for the m constraints -y_i <= 0. % zet = Lagrange multiplier for the single constraint -z <= 0. % s = Slack variables for the m general MMA constraints. % low = Column vector with the lower asymptotes, calculated and used % in the current MMA subproblem. % upp = Column vector with the upper asymptotes, calculated and used % in the current MMA subproblem. % %epsimin = sqrt(m+n)*10^(-9); epsimin = 10^(-7); raa0 = 0.00001; albefa = 0.1; asyinit = 0.5; %Sinit asyincr = 1.2; %Sfaster asydecr = 0.7; %Sslower move = 0.2; %move = 1.0; %one can try for different move value eeen = ones(n,1); eeem = ones(m,1); zeron = zeros(n,1); %%% these lines are required to reshape so that matrix dimensions %%% satisfied for solving four, nine and six node problems with MMA. %xval = reshape(xval,m*n,1); %xold1 = reshape(xold1,m*n,1); %xold2 = reshape(xold2,m*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); %10 lowmax = xval - 0.01*(xmax-xmin); uppmin = xval + 0.01*(xmax-xmin); uppmax = xval + 10*(xmax-xmin); %10 low = max(low,lowmin); low = min(low,lowmax); upp = min(upp,uppmax); upp = max(upp,uppmin); end % Calculation of the bounds alfa and beta : 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); % Calculations of p0, q0, P, Q and b. xmami = xmax-xmin; xmamieps = 0.00001*eeen; xmami = max(xmami,xmamieps); xmamiinv = eeen./xmami; ux1 = upp-xval; ux2 = ux1.*ux1; %(Uk-xk)^2 xl1 = xval-low; xl2 = xl1.*xl1; %(xk-Lk)^2 uxinv = eeen./ux1; %1./(Uk-xk) xlinv = eeen./xl1; %1./(xk-Lk) % p0 = zeron; q0 = zeron; p0 = max(df0dx,0); q0 = max(-df0dx,0); %%% these lines are required to reshape so that matrix dimensions %%% satisfied for solving four,nine and six node problems with MMA. %p0 = reshape(p0,m*n,1); %q0 = reshape(q0,m*n,1); %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; % %%% Solving the subproblem by a primal-dual Newton method [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);