miqp.m

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

   end			     
end

% checking whether the bounds 0,1 on the binary variables are already present in
% the problem constraints, if not, add them

aux1         = lb(vartype);
index1       = find(aux1<0);
aux1(index1) = 0;
lb(vartype)  = aux1;

aux2         = ub(vartype);
index2       = find(aux2>1);
aux2(index2) = 1;
ub(vartype)  = aux2;
                                              
% The variable STACK is used to store the relaxed QP problems, that are
% generated during the Branch and Bound algorithm. 
% It's global to allow the subroutines at the end of the m-file to access it.
% STACKSIZE denotes the number of subproblems on the stack.          
% For almost all braching 
% strategies STACK has indeed the purpose of a stack (last in, first out).
% For the breadth first strategy however it turns out, that a data structure
% of a queue (first in, first out) is taken on. 

global STACK
global STACKSIZE
global STACKCOST

% Initialization of STACK with the MIQP                                

STACKSIZE = 1;
STACK     = struct('H',H, 'f',f,   'A',A,   'b',b,   'Aeq',Aeq, 'beq',beq, ...
                   'e',0, 'lb',lb, 'ub',ub, 'x0',x0, 'vartype',vartype, ...
		   'ivalues',-ones(nivar,1), 'level',0, 'xxmax', xxmax );

% the parameters in the structure STACK are:
% H,f         : parameters of the cost function
% A,b,Aeq,beq : parameters of the constraints
% e           : parameter to determine the cost of the optimization problem
% lb,ub,x0,vartype : parameters of the problem as in the original problem
% ivalues : values of the integer variables (vector of length length(vartype) )
% level   : depth of the node in the tree
  
STACKCOST = 0;    % Array storing the cost of the father problem, ordered in
                  % decreasing fashion (STACKCOST(1)=largest value)

% ==============================================================================
% Main Loop
% ==============================================================================

while (STACKSIZE > 0)&(QPiter < maxqp)&(flag ~= -1)
  
   % Get the next subproblem from the STACK                            
   subprob = pop;
  
   % Solve the relaxed qp or lp
   if size(subprob.H,1)>0
   
      switch solver

         case 'lp'
	    % collects all constraints in one matrix and one vector such that
	    % the first constraints are equality constraints. Calling lp the
	    % first rows of the constraints are declared as equality constraints
          
	    swarn = warning;
	    warning off;
            [x,lambda,how] = lp(subprob.f, [subprob.Aeq; subprob.A], ...
	     [subprob.beq; subprob.b], subprob.lb, subprob.ub, subprob.x0, ...
	     size(subprob.Aeq,1));
	    warning(swarn); 

         case 'linprog'

	    swarn = warning;
	    warning off;
            [x, fval, exitflag, outpu, lambda] = linprog(subprob.f, ...
             subprob.A, subprob.b, subprob.Aeq, subprob.beq, subprob.lb,...
             subprob.ub, subprob.x0, optlinprog);        
	    warning(swarn); 
            if exitflag > 0
               how = 'ok'; 
            elseif exitflag == 0
               warning('maximum number of iterations occurred')
               how = 'ok';  
            else
               how = 'infeasible';
	       if verbose >= 2
                  warning('no distinction whether unbounded or infeasible') 
	       end	  
            end 

         case 'lpnag'
	    if ~isempty(Aeq)
	       warning('Equality constraints are not explicitly supported')
	       warning('Will switch to Aeq x <= beq, -Aeq x <= -beq')
	    end
	    blpnag = [ subprob.b; subprob.beq; -subprob.beq ];
	    Alpnag = [ subprob.A; subprob.Aeq; -subprob.Aeq ];
            ifail  = 1;
            bl     = [subprob.lb; -inftol*ones(length(blpnag),1)];  
            bu     = [subprob.ub;  blpnag]; 
	               
	    swarn = warning;
 	    warning off;
            [x, istate, objlp, clamda, ifail] = e04mbf(bl, bu, subprob.x0, ...
             subprob.f, Alpnag, verbose-1, maxQPiter, ifail);
	    warning(swarn); 
	     
            switch ifail
               case {-1,0}
                  how = 'ok';
               case {2}
                  how = 'unbounded';
               case {3,4}
                  % might also be considered as infeasible 
                  warning('QP is cycling or too few iterations')
                  how = 'ok';
               case {1}
                  how = 'infeasible';
               otherwise
                  error('other error code in ifail from e04mbf')
            end	     

         case 'qp'
	    if isfield(Options,'postol')
	       if verbose >= 1
	          if cond(H) <= postol
		     warning('H is close to singularity')
		  end   
	       end
	    end
	    % collects all constraints in one matrix and one vector such that
	    % the first constraints are equality constraints. Calling lp the
	    % first rows of the constraints are declared as equality constraints
	    
	    swarn = warning;
 	    warning off;
	    [x,lambda,how] = qp(subprob.H, subprob.f, ...
	     [subprob.Aeq; subprob.A], [subprob.beq; subprob.b], ...
             subprob.lb, subprob.ub, subprob.x0, size(subprob.Aeq,1));
	    warning(swarn); 

         case 'quadprog'
	    if isfield(Options,'postol')
	       if verbose >= 1
	          if cond(H) <= postol
		     warning('H is close to singularity')
		  end   
	       end
	    end

	    swarn = warning;
 	    warning off;
            [x,fval,exitflag,outpu,lambda] = quadprog(subprob.H, subprob.f, ...
             subprob.A, subprob.b, subprob.Aeq, subprob.beq, subprob.lb, ...
             subprob.ub, subprob.x0, optquadprog);        
	    warning(swarn); 

            if exitflag > 0
               how = 'ok'; 
            elseif exitflag == 0
               warning('maximum number of iterations occurred')
               how = 'ok';  
            else
               how = 'infeasible';
	       if verbose >= 2
                  warning('no distinction whether unbounded or infeasible')
	       end	  		  
            end

         case 'qpnag'
	    if isfield(Options,'postol')
	       if verbose >= 1
	          if cond(H) <= postol
		     warning('H is close to singularity')
		  end   
	       end
	    end
	    if ~isempty(Aeq)
	       warning('Equality constraints are not explicitly supported')
	       warning('Will switch to Aeq x <= beq, -Aeq x <= -beq')
            end
	    bqpnag = [ subprob.b; subprob.beq; -subprob.beq ];
	    Aqpnag = [ subprob.A; subprob.Aeq; -subprob.Aeq ];	  	    
            bl     = [subprob.lb; -inftol*ones(length(bqpnag),1)];  
            bu     = [subprob.ub;  bqpnag];            
            istate = zeros(length(bu),1);
            featol = wu*ones(length(bu),1);
            ifail  = 1;
	    
	    swarn = warning;
 	    warning off;
            [x, iter, obj, clamda, istate, ifail] = ...
             e04naf(bl, bu, 'qphess', subprob.x0, subprob.f, Aqpnag, ...
             subprob.H, lpf, cold,istate, featol, verbose-1, maxQPiter, ...
	     inftol, orthog, ifail); 
	    warning(swarn); 
            
	    switch ifail
               case {0,1,3}
                  how = 'ok';
               case 2
                  how = 'unbounded';
               case {4,5}
                  % might also be considered as infeasible 
                  warning('QP is cycling or too few iterations')
                  how = 'ok';
               case {6,7,8}
                  how = 'infeasible';
               otherwise
                  error('other error code in ifail from e04naf')
            end
            
         case 'qp_dantz'
         
            if ~isempty(subprob.Aeq)
               error('this solver does not support equality constraints')
            end   
            
            [x, how] = qp_dantz( subprob.H,   subprob.f, ...
                [subprob.A; eye(length(subprob.f));-eye(length(subprob.f))], ...
                [subprob.b; subprob.ub(:);         -subprob.lb(:)], ...
                 subprob.xxmax(:));
            
            if strcmp(how,'feasible')
               how = 'ok';
            end             

         otherwise
            error('unknown solver')

      end
      QPiter = QPiter + 1;
   else
      x=[];
      if all(subprob.b>=0),
         how = 'ok';
      else
         how = 'infeasible';
      end
   end
   
   if strcmp(how,'unbounded')
   % If the relaxed problem is unbounded, so is the original problem 
      xmin = [];
      fmin = -inf;
      flag = -1;      
      warning('unbounded cost function')
      return
   elseif strcmp(how,'infeasible')
      % subproblem fathomed                                       
   else   
      % subproblem feasible                                      
     
      if ~isempty(x), 
         zpi=.5*x'*subprob.H*x+x'*subprob.f+subprob.e;
      else
         zpi=subprob.e;
      end
      if flag~= 1   % if no integer feasible solution has been found yet,
         flag = 5;  % the problem is now at least feasible
      end
     
      % Check if value function is better than the value so far    
      if zpi<=zstar
                           
         xi               = x(subprob.vartype); % integer variables
         xc               = x;
         xc(subprob.vartype) = [];              % continuous variables
              
         % Test whether the relaxed integer variables are feasible,  
         % i.e. whether they are integral. Note that this condition  
         % is always satisfied, if there are no free integer variables,     
         % i.e if vartype is empty  
	  
         if norm( round(xi) - xi,inf) < integtol
            % subproblem solved                                   
            % update the value of the cost function              
            zstar             = zpi;
            ifree             = find(subprob.ivalues==-1);
            iset              = find(subprob.ivalues>-1);
            absi              = 1:nx;
            absi(cont)        = [];
            if rounding
               xstar(absi(ifree))= round(xi);
               xstar(absi(iset)) = round(subprob.ivalues(iset));
            else
               xstar(absi(ifree))= xi;
               xstar(absi(iset)) = subprob.ivalues(iset);
            end 
            xstar(cont)       = xc;
            flag              = 1;
	    optQP             = QPiter;
         else
            % separate subproblem, if there are still free integer   
            % variables. Note that no further test is required,      
            % whether there are still integer variables, since ixrel 
            % is nonempty                                            
            % branchvar denotes the position of the branching        
            %           integer variable within the set of integer 
            %           variables in the present subproblem          
            branchvar            = decision(x, subprob.vartype, branchrule);
            [p0,p1,zeroOK,oneOK] = separate(subprob, branchvar);
            switch method
               case 'depth'
                  cost = 1/(subprob.level+1); 
               case 'breadth'
                  cost = subprob.level+1;
               case 'best'
                  cost = zpi; % Best-first. This tends to go breadth-first
               case 'bestdepth' 
                  cost = zpi/(subprob.level+1); % This privilegiates deep nodes
            end
            if order == 0
               if oneOK
                  push(p1,cost);
               end
               if zeroOK   
                  push(p0,cost);
               end
            else
               if zeroOK
                  push(p0,cost);
               end
               if oneOK
                  push(p1,cost);
               end   
            end  %if order ...
         end  %if norm ... 
      end  %if zpi<=zpstar ...
   end  %if strcmp ...

   if (QPiter>=maxqp)
      if (flag == 1) | (flag == 5)
         flag = flag+10; % update flag in case of limit of QPs reached
      end
   end   
   
   % Display present status of the MIQP                               
   if verbose >= 1
      disp('QPiter = ') , disp(QPiter)
      disp('zstar  = ') , disp(zstar)
      disp('xstar  = ') , disp(xstar')
   end