lipsol.m

数学建模的源代码

   
   [trerror,rrb,rrf,rru,rdgap,fval,dualval] = ...
      errornobj(b,rb,bnrm,f,rf,fnrm,ub,ru,unrm,dgap,Ubounds_exist,x,y,w);
   
   iter = 0;
   converged = 0;
   backed = 0;
   if (~Ubounds_exist)
      sn1 = [];
   end
   
   if (diagnostic_level >= 2)
      if (Ubounds_exist)
         fprintf(['\n  Residuals:   Primal     Dual     Upper    Duality     Total']);
         fprintf(['\n               Infeas    Infeas    Bounds     Gap        Rel\n']);
         fprintf(['               A*x-b   A''*y+z-w-f {x}+s-ub  x''*z+s''*w   Error\n']);
         fprintf('  -------------------------------------------------------------\n');
      else
         fprintf(['\n  Residuals:   Primal     Dual     Duality    Total']);
         fprintf(['\n               Infeas    Infeas      Gap       Rel\n']);
         fprintf(['               A*x-b    A''*y+z-f    x''*z      Error\n']);
         fprintf('  ---------------------------------------------------\n');
      end
   end
   
   if (showstat)
      legstr = {'duality gap','primal infeas','dual infeas'};
      if (Ubounds_exist)
         legstr{1,4} = 'upper bounds';
      end
   end
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Begin main loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   while (iter <= maxiter)
      
      if (diagnostic_level >= 2)
         fprintf('  Iter %4i:  ', iter);
         if (Ubounds_exist)
            fprintf('%8.2e %8.2e %8.2e %8.2e %8.2e\n',rb,rf,ru,dgap,trerror);
         else
            fprintf('%8.2e %8.2e %8.2e %8.2e\n',rb,rf,dgap,trerror);
         end
      end
      if (showstat)
         if (iter == 0)
            lipsolfig = figure;
            semilogy(0,dgap,'c-',0,rb,'b:',0,rf,'g-.');
            dg_prev = dgap;
            rb_prev = rb;
            rf_prev = rf;
            hold on
            if (Ubounds_exist)
               semilogy(ru,'r--');
               ru_prev = ru;
            end
            drawnow
            title(['linprog']); 
            xlabel('Iteration Number');
            ylabel('Residuals');
         else % (iter > 0)
            figure(lipsolfig)
            semilogy([iter-1 iter],[dg_prev dgap],'c-', ...
               [iter-1 iter],[rb_prev rb],'b:', ...
               [iter-1 iter],[rf_prev rf],'g-.');
            dg_prev = dgap;
            rb_prev = rb;
            rf_prev = rf;
            if (Ubounds_exist)
               semilogy([iter-1 iter],[ru_prev ru],'r--');
               ru_prev = ru;
            end
            set(gca,'XLim',[0 iter+1])
            legend(legstr{:})
            drawnow
         end
      end
      
      if (iter > 0)
         stop = 0;
         converged = 0;
         if (diagnostic_level >= 1)
            message = '';
         end
         trerror = Hist(1,iter);
         if (trerror < tol)
            stop = 1;
            converged = 1;
            if (diagnostic_level >= 1)
               message = 'Optimization terminated successfully.';
            end
         else
            small = 1.e-5;
            if (iter > 2)
               blowup = (Hist(1:5,iter) > (max(tol,min(Hist(1:5,:)')')/small));
               if (any(blowup))
                  stop = 1;
                  converged = 0;
                  if (diagnostic_level >= 1)
                     message ...
                     = sprintf(['One or more of the residuals, duality gap, or total relative error\n', ...
                        '         has grown 100000 times greater than its minimum value so far:\n']);
                     message = [message detectinf(tol,Hist)];
                  end
               end
               nstall = 3;
               if (iter > nstall)
                  latest = iter-nstall:iter;
                  h = Hist(1:5,latest);
                  hmean = mean(h');
                  for i = 1:5
                     stall(i) = ((hmean(i) > small) & ...
                        (all(abs(h(i,:)-hmean(i)) < small*hmean(i))));
                  end
                  if (any(stall))
                     stop = 1;
                     converged = 0;
                     if (diagnostic_level >= 1)
                        message ...
                        = sprintf(['One or more of the residuals, duality gap, or total relative error\n', ...
                              '         has stalled:\n']);
                        message = [message detectinf(tol,Hist)];
                     end
                  end
               end
            end
         end
         if (stop)
            break % out of while (iter <= maxiter) loop
         end
      end % if (iter > 0)
      
      xn1 = reciprocal(x);
      if (Ubounds_exist)
         sn1 = reciprocal(s);
         vmid = reciprocal(z.*xn1 + w.*sn1);
      else
         vmid = reciprocal(z.*xn1);   
      end
      
      vmid = full(min(1e+15,vmid));
      
      [P,U] = getpu(A,vmid,Dense_cols_exist,idense,ispars);   
      
      [dx,dy,dz,ds,dw,Sherman_OK,Mdense,Rinf,cgiter] = ...
         direction(A,P,U,Rb,Rf,Ru,Rxz,Rsw,vmid,xn1,sn1,z,w,0,1, ...
         Sherman_OK,Dense_cols_exist,Ubounds_exist,Mdense,perm,Rinf,Hist,cgiter);
      
      [ap,ad] = ratiotest(dx,dz,ds,dw,Ubounds_exist,x,z,s,w);
      
      if ((tau0*ap < 1) | (tau0*ad < 1))
         
         newdgap = (x + min(1,ap)*dx)'*(z + min(1,ad)*dz);
         if (Ubounds_exist) 
            newdgap = newdgap + (s + min(1,ap)*ds)'*(w + min(1,ad)*dw); 
         end
         sigmak = (newdgap/dgap)^2;
         sigmin = 0;
         sigmax = .208; % Zhang's choice
         p = ceil(log10(trerror)); 
         if ((p < -1) & (dgap < 1.e+3))
            sigmax = 10^(p+1);
         end
         sigmak = max(sigmin, min(sigmax, sigmak));
         mu = sigmak*dgap/nt;      
         Rxz = dx.*dz;
         Rsw = ds.*dw;
         
         [dx2,dy2,dz2,ds2,dw2,Sherman_OK,Mdense,Rinf,cgiter] = ...
            direction(A,P,U,sparse(m,1),sparse(n,1),sparse(n,1),...
            Rxz,Rsw,vmid,xn1,sn1,z,w,mu,2,Sherman_OK, ...
            Dense_cols_exist,Ubounds_exist,Mdense,perm,Rinf,Hist,cgiter);
         
         dx = dx + dx2;
         dy = dy + dy2;
         dz = dz + dz2;
         ds = ds + ds2;
         dw = dw + dw2;
         
         [ap,ad] = ratiotest(dx,dz,ds,dw,Ubounds_exist,x,z,s,w);
         
      end
      
      tau = tau0; 
      if (~backed)
         tau = .9 + 0.1*tau0;
      end
      k = ceil(log10(trerror));
      if (k <= -5)
         tau = max(tau,1-10^k);
      end
      ap2 = min(tau*ap,1);
      ad2 = min(tau*ad,1);
      xc = x;
      yc = y;
      zc = z;
      sc = s;
      wc = w;
      step = [1 .9975 .95 .90 .75 .50];
      for k = 1:length(step)
         x = xc + step(k)*ap2*dx;
         y = yc + step(k)*ad2*dy;
         z = zc + step(k)*ad2*dz;
         s = sc + step(k)*ap2*ds;
         w = wc + step(k)*ad2*dw;
         [Rxz,Rsw,dgap] = complementy(x,z,s,w,Ubounds_exist);
         phi = nt*full(min([Rxz;Rsw(find(Rsw))]))/dgap;
         if ((max(ap2,ad2) == 1) | (phi >= phi0))
            break
         end
      end
      phi0 = min(phi0, phi);
      if (k > 1)
         backed = 1;
      end
      
      [Rb,Rf,rb,rf,ru,Ru] = feasibility(A,b,f,ub,Ubounds_exist,x,y,z,s,w);
      if any(isnan([rb rf ru]))
         xsol = [];
         fval = [];
         lambda.ineqlin = [];
         lambda.eqlin = [];
         lambda.upper = [];
         lambda.lower = [];
         exitflag = -1;
         output.iterations = iter;
         output.cgiter = cgiter;
         return
      end
      
      [trerror,rrb,rrf,rru,rdgap,fval,dualval] = ...
         errornobj(b,rb,bnrm,f,rf,fnrm,ub,ru,unrm,dgap,Ubounds_exist,x,y,w);
      
      Hist = [Hist [trerror rrb rrf rru rdgap fval dualval]'];
      iter = iter + 1;
      
   end % while (iter <= maxiter)
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End main loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
else % All variables are fixed, singelton, etc - no more variables left to solve for.
   iter = 0;
   cgiter = 0;
   message = '';
   converged = 1;
   x = zeros(0,1);
   w = [];
   y = [];
   s = [];
   z = [];
end

xsol = x;

if (col_scaled)
   xsol = xsol .* colscl;
   if (diagnostic_level >= 4)
      fprintf(['  Scaling solution back to original problem.\n']);
   end
end

if (Lbounds_non0)
   xsol = xsol + lb;
   if (diagnostic_level >= 4)
      fprintf(['  Shifting solution back away from zero to reflect' ...
            ' original lower bounds.\n']);
   end
end

if (nbdslack < n)
   % Remove slack solution variables
   xsol = xsol(1:nbdslack);
   if (diagnostic_level >= 4)
      fprintf(['  Eliminating %d slack variables (bounds) from solution.\n'], ...
         n-nbdslack);
   end
end

if (Sgtons_exist)
   xtmp(insolved,1) = xsol;
   xtmp(isolved,1) = xsolved;
   xsol = xtmp; 
   n = length(xsol);
   if (diagnostic_level >= 4)
      fprintf(['  Reuniting %d singleton variables with solution.\n'], ...
         length(xsolved));
   end
end

if (Zrcols_exist)
   xtmp(inzcol,1) = xsol;
   xtmp(izrcol,1) = xzrcol;
   xsol = xtmp;
   n = length(xsol);
   if (diagnostic_level >= 4)
      fprintf(['  Reuniting %d obvious zero variables with solution.\n'], ...
         length(xzrcol));
   end
end

if (Row_deleted)
   if (norm(Adel*xsol-bdel)/max(1,norm(bdel)) > tol)
      converged = 0;
      if (diagnostic_level >= 1)
         message = sprintf(['the primal is infeasible.\n' ...
         '         The equality constraints are dependent but not consistent.']);
      end
   else
      if (diagnostic_level >= 4)
         fprintf(['  Solution satisfies %d dependent constraints.\n'], ...
            size(Adel,1));
      end
   end
end

if (Fixed_exist)
   xtmp(infx,1) = xsol;
   xtmp(ifix,1) = xfix;
   xsol = xtmp;
   n = length(xsol);
   if (diagnostic_level >= 4)
      fprintf(['  Reuniting %d fixed variables with solution.\n'], ...
         length(ifix));
   end
end

if (nslack < n)
   % Remove slack solution variables
   xsol = xsol(1:nslack);
   if (diagnostic_level >= 4)
      fprintf(['  Eliminating %d slack variables from solution.\n'], ...
         n-nslack);
   end
end

if (nlbinf > 0) % if there were unbounded below variables split up
   xsol(lbinf) = xsol(lbinf) - xsol(n_orig+1:end);
   xsol = xsol(1:n_orig);
   if (diagnostic_level >= 4)
      fprintf(['  Subtracting the strictly positive differences of %d' ...
            ' unbounded below variables.\n'], ...
         nlbinf);
   end
end

fval = full(f_orig' * xsol);

if computeLambda
   % Calculate lagrange multipliers
   lbindexfinal = lbindex(lbindex <= lbend); 
   ubindexfinal = ubindex(ubindex <= ubend);
   zindexfinal = find(zindex);
   windexfinal = find(windex);
   nnzlbindexfinal = nnz(lbindexfinal);
   nnzubindexfinal = nnz(ubindexfinal);
   if boundsInA
      lambda.upper(ubindexfinal) = w(windexfinal(1:nnzubindexfinal)) - y(1:nnzubindexfinal);
      lambda.lower(lbindexfinal) = z(zindexfinal(1:nnzlbindexfinal)) - ...
         y(nnzubindexfinal+1:nnzlbindexfinal+nnzubindexfinal);  
      if Fixed_exist 
         lambda.lower(fixedlower) = f_before_deletes(fixedlower);
         lambda.upper(fixedupper) = -f_before_deletes(fixedupper);
      end
      if Sgtons_exist
         % sgrows: what eqlin lambdas we are solving for (row in A)
         % sgcols: what column in A that singleton is in
         result = -f_orig + lambda.lower - lambda.upper;