lipsol.m

线性规划算法

   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])) 
          if (diagnostic_level >= 1) 
              fprintf(['Exiting: cannot converge because the primal residual, dual residual,  \n' ... 
                      '         or upper-bound feasibility is NaN.\n']); 
          end 
          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 = 'Solution determined by constraints.';     
   converged = 1; 
   x = zeros(0,1); 
   y = zeros(0,1); 
   z = zeros(0,1); 
   w = zeros(0,1); 
   s = []; 
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; 
         lambda.eqlin(sgrows) = result(sgcols); 
      end 
   else 
      iAineqindex = Aindex( find (Aindex >= Aineqstart & Aindex <= Aineqend) ); 
      % Subtract off extraAineqvars and numRowsAineq even if those were