%-------------------------------------------------------
%    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);