panmip.m

由matlab开发的hybrid系统的描述语言

      G           = H;
      bl          = [vlb; -huge*ones(nlin,1)];
      bu          = [vub;  d];
      nivar       = length(ivar);
      ivar        = ivar(:);

      [xopt, fopt, flag] = miqp_d(a_dense, nvar, nlin, x0, problemtype,...
                G, bl, bu, nivar, ivar, priority, iprint, epsilon); 
                      
      flag = flag(length(flag));
      
   case {'sunbaronQP','sunbaronLP'}
   %-------------------------------      
      dircnt       = dir;
      optionsthere = 0;
      ii = 1;
      if strcmp(Options.solver,'sunbaronQP')
         while ((ii <= size(dircnt,1)) & ~optionsthere)
            if strcmp('optionsQP',dircnt(ii).name)
               optionsthere = 1;
            end
            ii = ii+1;
         end
         if ~optionsthere
            dircnt       =   dir(BaronDirectory);
            optionsthere = 0;
            ii = 1;    
            while ((ii <= size(dircnt,1)) & ~optionsthere)
               if strcmp('optionsQP',dircnt(ii).name)
                  optionsthere = 1;
               end
               ii = ii+1;
            end
            if ~optionsthere          
               error('Baron requires file optionsQP in the working directory')
            else
               loactionsOQP = strcat(BaronDirectory,'/optionsQP');
               statuscp     = copyfile(loactionsOQP,'options'); 
               if ~statuscp
                  error('failed to copy the file optionsQP')
               end
            end   
         else
            copyfile('optionsQP','options')
         end   
      else
         while ((ii <= size(dircnt,1)) & ~optionsthere)
            if strcmp('optionsLP',dircnt(ii).name)
               optionsthere = 1;
            end
            ii = ii+1;
         end
         if ~optionsthere
            dircnt       =   dir(BaronDirectory);
            optionsthere = 0;    
            ii = 1;    
            while ((ii <= size(dircnt,1)) & ~optionsthere)
               if strcmp('optionsLP',dircnt(ii).name)
                  optionsthere = 1;
               end
               ii = ii+1;
            end
            if ~optionsthere          
               error('Baron requires file optionsLP in the working directory')
            else 
               loactionsOLP = strcat(BaronDirectory,'/optionsLP');
               statuscp     = copyfile(loactionsOLP,'options');
               if ~statuscp
                  error('failed to copy the file optionsLP')
               end
            end   
         else
            copyfile('optionsLP','options')
         end
      end               
      
      [A, b] = ctype2matCompact(C, d, ctype, rangevar);
      ivar   = vartype2ivar(vartype, n, 'printerror');
      [lb,ub]= add01(ivar, lb, ub);

      if strcmp(Options.solver,'sunbaronQP')
         if ~isfield(Options,'sunbaronQP')
            sunbaroptions = [];
            infty         = 1e8;
         else
            sunbaroptions = Options.sunbaronQP;
	    if isfield(Options.sunbaronQP,'infinity')
               infty = Options.sunbaronQP.infinity;
            else
               infty = 1e8;
            end   
         end   
         baronlpflag = 0;
      else
         if ~isfield(Options,'sunbaronLP')
            sunbaroptions = [];
            infty         = 1e8;
         else
            sunbaroptions = Options.sunbaronLP;
	    if isfield(Options.sunbaronLP,'infinity')
               infty = Options.sunbaronLP.infinity;
            else
               infty = 1e8;
            end   
         end         
         baronlpflag = 1;
      end

      if strcmp(Options.solver,'sunbaronLP')
         if ~isfield(Options,'sunbaronLP')
            suppresserr = 0;	 	     
         elseif isfield(Options.sunbaronLP,'suppresserr') 
            suppresserr = Options.sunbaronLP.suppresserr;	 
         else
            suppresserr = 0;	 
	 end   
      else
         suppresserr = 0;
      end	 	
      	
      % remove infinity terms in box constraints
      for i=1:length(lb)
         if lb(i) == -inf
            lb(i) = -infty;
         end
         if ub(i) == inf 
            ub(i) = infty;
         end
      end             
      
      if ( (norm(H,inf) > 1e-4) & strcmp(Options.solver,'sunbaronLP') )
         if ~suppresserr
            disp('This problem has a nonzero H and cannot be solved with MILP')
            error('Use the solver sunbaronQP instead of sunbaronLP. Thanks!')
	 else
	    H = 0*H; 
	 end      
      end
            
      [xopt, fopt, flag, qpc ] = ...
       baronsolver( H, f, A, b, lb, ub, ivar, sunbaroptions, baronlpflag, ...
                    BaronDirectory);
 
      delete options % clean up directory:optionsLP or optionsQP are still there
      
      if (flag == 7) | (flag == 0)
         xopt = x_opt_init;
         fopt = f_opt_init;
      end                         

       
   case {'mps'}
   %-----------  
      if ~isfield(Options,'mps')
         Options.mps = [];
      end
                  
      Al = [];  % Al x <= bl 
      bl = [];
      Ag = [];  % Ag x >= bg
      bg = [];
      Ae = [];  % Ae x == be
      be = [];
      Ar = [];  % ranged constraints
      br = [];
      
      if isfield(Options.mps,'extractbox')
         extractbox = Options.mps.extractbox;
      else
         extractbox = 0;         
      end
      if isfield(Options.mps,'verbose')
         verbose = Options.mps.verbose;
      else
         verbose = 0;         
      end
      if isfield(Options.mps,'outname')
         prname = Options.mps.outname;
      else
         prname = 'mpsfile.mps';
      end  
      if isfield(Options.mps,'infinity')
         infty  = Options.mps.infinity;
      else
         infty  = 1e8;
      end   
      if isfield(Options.mps,'format')
         format  = Options.mps.format;
      else
         format  = 'cplex';
      end              
      if isfield(Options.mps,'QPflag')
         QPflag  = Options.mps.QPflag;
      else
         QPflag  = 0;
      end 
      if isfield(Options.mps,'LPflag')
         LPflag  = Options.mps.LPflag;
      else
         LPflag  = 0;
      end
      if isfield(Options.mps,'MILPflag')
         MILPflag  = Options.mps.MILPflag;
      else
         MILPflag  = 0;
      end
      if isfield(Options.mps,'suppresserr')
         suppresserr  = Options.mps.suppresserr;
      else
         suppresserr  = 0;
      end 
                      
      lower    = lb;
      upper    = ub;
      
      if extractbox
         todelete = [];       
         i        = 1;
     
         % extract contraints from the general constraint matrix, which are
         % actually box bounds
      
         while i <= size(C,1)
            [xi,yi,vi] = find(C(i,:));
            if length(xi)==1  % there is exactly one nonzero element in i-th row
               if strcmp(ctype(i),'L')
                  if vi > 0                    % vi is nonzero by construction
                     if upper(yi)>d(i,:)/vi    % the bound in C is tighter than
                        upper(yi) = d(i,:)/vi; % the previous box bound
                     end
                  else
                     if lower(yi) < d(i,:)/vi
                        lower(yi) = d(i,:)/vi;
                     end
                  end
                  todelete = [todelete; i];  % collect indices to delete
               elseif strcmp(ctype(i),'G')
                  if vi > 0
                     if lower(yi) < d(i,:)/vi 
                        lower(yi) = d(i,:)/vi;
                     end
                  else      
                     if upper(yi) > d(i,:)/vi
                        upper(yi) = d(i,:)/vi;
                     end
                  end 
                  todelete = [todelete; i];  % collect indices to delete
               end                         
            end                      
            i = i+1;
         end
      
         C(todelete,:)   = []; 
         d(todelete,:)   = [];
         ctype(todelete) = [];
      end
      
      ivar             = vartype2ivar(vartype, n, 'noerror');
      [lower, upper]   = add01(ivar, lower, upper);
      
      [Al,bl, Ag,bg, Aeq,beq, Ar,br, r2] = ctype2allmat(C,d,ctype,rangevar);

      ivartrue = vartype2integer(vartype,n);
            
      writemps(prname, H,f, Al,bl, Ag,bg, Aeq,beq, Ar,br, r2, ... 
               lower,upper, ivar,ivartrue, infty, format,verbose, ...
               QPflag, MILPflag, LPflag, suppresserr);
      
      xopt = x_opt_init;
      fopt = f_opt_init;
      flag = 2;
      
      if verbose >=1
         disp('mps file generated, no values assigned to output parameters')
      end         

   case {'cplex'}
   % ------------
      if ~isfield(Options,'cplex')
         Options.cplex = [];
      end
   
      if ~isfield(Options.cplex,'verbose')
         verbose = 0;
      else
         verbose = Options.cplex.verbose;   
      end  
          
      if ~isfield(Options.cplex,'savecplex')
         savecplex = 0;
      else
         savecplex = Options.cplex.savecplex;   
      end      
      
      if ~isempty(x0)
         x0i = 1:length(f);
      else
         x0i = [];
      end    
      
      ivar   = vartype2ivar(vartype,n,'noerror');
      [lb,ub]= add01(ivar,lb,ub);   
      
      [xopt, fopt, lpsolved, status, slack]= milp_cplex(1, f, C, d, ctype(:),...
      lb(:), ub(:), rangevar(:), vartype(:), x0i(:), x0(:), verbose, savecplex);
       
      Extendedflag.lpsolved = lpsolved;
      Extendedflag.status   = status;
      Extendedflag.slack    = slack;
                                         
      switch status
         case {1}   % feasible
            flag = 1;
         case {0,2} % infeasible