function [x,fval,exitflag,output,lambda,nlp,pos]=milprog(f,n,m,A,B,Aeq,Beq,lb,ub,x0,options)
// 
// //MILPROG     Linear programming by exaustive search.                 
// //   X=MILPROG(f,n,m,A,b) solves the mixed integer linear programming problem:
// //        
// //            min f'*X    subject to:   A*X <= b 
// //            x,y
// //            x real
// //            y integer
// //            X=[x;y]
// //            n=number of real variables
// //            m=number of integer variables
// //   X=MILPPROG(f,n,m,A,b,Aeq,beq) solves the problem above while additionally
// //   satisfying the equality constraints Aeq*[x;y] = beq.
// //   
// //   X=MILPROG(f,n,m,A,b,Aeq,beq,LB,UB) defines a set of lower and upper
// //   bounds on the design variables, real x, so that the solution is in
// //   the range LB <= x <= UB.  Use empty matrices for LB and UB
// //   if no bounds exist. Set LB(i) = -Inf if x(i) is unbounded below; 
// //   set UB(i) = Inf if x(i) is unbounded above.
// //
// //   X=MILPROG(f,n,m,A,b,Aeq,beq,LB,UB,x0) sets the starting point of real
// //   variables to x0.  This option is only available with the active-set algorithm.  The default
// //   interior point algorithm will ignore any non-empty starting point.
// //
// //   X=MILPROG(f,n,m,A,b,Aeq,Beq,LB,UB,x0,OPTIONS) minimizes with the default 
// //   optimization parameters replaced by values in the structure OPTIONS, an 
// //   argument created with the OPTIMSET function.  See OPTIMSET for details.  
// //   Use options are Display, Diagnostics, TolFun, LargeScale, MaxIter. 
// //   Currently, only 'final' and 'off' are valid values for the parameter 
// //   Display when LargeScale is 'off' ('iter' is valid when LargeScale is 'on').
// //
// //   [X,FVAL]=MILPROG(f,n,m,A,b) returns the value of the objective function at X:
// //   FVAL = f'*X.
// //
// //   [X,FVAL,EXITFLAG] = MILPROG(f,A,b) returns EXITFLAG that 
// //   describes the exit condition of LINPROG.
// //   If EXITFLAG is:
// //      > 0 then LINPROG converged with a solution X.
// //      0   then LINPROG reached the maximum number of iterations without converging.
// //      < 0 then the problem was infeasible or LINPROG failed.
// //
// //   [X,FVAL,EXITFLAG,OUTPUT] = MILPROG(f,n,m,A,b) returns a structure
// //   OUTPUT with the number of iterations taken in OUTPUT.iterations, the type
// //   of algorithm used in OUTPUT.algorithm, the number of conjugate gradient
// //   iterations (if used) in OUTPUT.cgiterations.
// //
// //   [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,NLP,POS]=MILPROG(f,n,m,A,b) returns the set of 
// //   Lagrangian multipliers LAMBDA, at the solution: LAMBDA.ineqlin for the 
// //   linear inequalities A, LAMBDA.eqlin for the linear equalities Aeq, 
// //   LAMBDA.lower for LB, and LAMBDA.upper for UB.NLP returns the number of
// //   LINPTOG evaluations.POS returns the position of the the solution, use to
// //   analize the results.
// //   
// //   NOTE: the LargeScale (the default) version of MILPROG uses LINPROG that use a primal-dual
// //         method. Both the primal problem and the dual problem must be feasible 
// //         for convergence. Infeasibility messages of either the primal or dual, 
// //         or both, are given as appropriate.  The primal problem in standard 
// //         form is 
// //              min f'*x such that A*x = b, x >= 0.
// //         The dual problem is
// //              max b'*y such that A'*y + s = f, s >= 0.
// // Authors: Argimiro R. Secchi &
// //          Marcelo Escobar
// //          Chemical Engineering Department
// //          Applied Math Department
// //          Universidade Federal do Rio Grande do Sul, Brasil
// //          arge@enq.ufrgs.br | escobar@enq.ufrgs.br
// //
// // version 1.0: 15/agosto/2004
// 
// If just 'defaults' passed in, return the default options in X

// If just 'defaults' passed in, return the default options in X

defaultopt = struct('Display','final',...
   'TolFun',1e-8,'Diagnostics','off',...
   'LargeScale','on','MaxIter',85);

if nargin==1 & nargout <= 1 & isequal(f,'defaults')
   x = defaultopt;
   return
end

// Handle missing arguments

 if nargin < 11, options = [];
   if nargin < 10, x0 = []; 
      if nargin < 9, ub = []; 
         if nargin < 8, lb = []; 
            if nargin < 7, Beq = [];
               if nargin < 6, Aeq = [];
               end, end, end, end, end, end
if nargout > 4
   computeLambda = 1;
else 
   computeLambda = 0;
end

//M-Matriz de Possibilidades:
nv=length(f);
ne=size(Aeq);
ni=size(A);
ny=m;
N=[];
M=zeros(2^ny,ny);

[x0,lb,ub,msg] = checkbounds(x0,lb,ub,n);
if ~isempty(msg)
   exitflag = -1;
   output = []; x=x0; fval = []; lambda = []; nlp = 0; pos = [];
   if verbosity > 0
      disp(msg)
   end
   return
end

for i=1:ny
   eq=[zeros(2^(ny-i),1);ones(2^(ny-i),1)];
   cp=length(eq);
      if cp==2^ny
          N=eq;
      else N=[];
          while cp~=2^ny
             N=[N;eq];
             cp=length(N);
             end
      end
 M(:,i)=N;
 end
 tol=1e-6;
 if nv==m // Apenas variaveis Binarias
     fv=[];
     pos=[];
    if ne(1)==0 //Apenas Equacoes de Desigualdade
         for i=1:2^ny
         flagv(i)=sum((A*M(i,:)'-B)>=0);   //flagv=1 implica em violacao de restricao
         if flagv(i)==0
             fv=[fv;f*M(i,:)'];
             pos=[pos;i];
         end
         end
     fval=min(fv);
     x=M(pos(find(fv==min(fv))),:)';
     pos=pos(find(fv==min(fv)));
    else //Equacoes de Igauladdae e Desigualdade
         for i=1:2^ny
         flagv(i)=sum(abs(Aeq*M(i,:)'-Beq)>=tol)+sum((A*M(i,:)'-B)>=0);  //flagv=1 implica em violacao de restricao
         if flagv(i)==0
             fv=[fv;f*M(i,:)'];
             pos=[pos;i];
         end
         end
     fval=min(fv);
     x=M(pos(find(fv==min(fv))),:)';
     pos=pos(find(fv==min(fv)));
    end
     exitflag=1;
     output=[];
     lambda=[];
     nlp=0;
else //Variaveis Binarias e Reais
    if ne(1)==0 //Apenas Equacoes de Desigualdade
        c=f(1:nv-m);
        Am=A(:,1:nv-m);
        lb=lb(1:nv-m);
        ub=ub(1:nv-m);
        op=options;
        fv=[];pos=[];Mxy=[];nlp=0;
        for i=1:2^ny
            Bm=B-f(nv-m+1:end)*M(i,:)';
            [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=linprog(c,Am,Bm,Aeq,Beq,lb,ub,x0,op);
            exitflag{i}=EXITFLAG;
            output{i}=OUTPUT;
            lambda{i}=LAMBDA;
            nlp=nlp+(EXITFLAG>0);
            xmy=[X;M(i,:)'];
            Mxy=[Mxy;xmy'];
            FVALy(i)=FVAL+f(nv-m+1:end)*M(i,:)';
            flagv(i)=sum((A*xmy-B)>=0);
            if flagv(i)==0
               fv=[fv;f*xmy];
               pos=[pos;i];
            end
        end
        fval=min(fv);
        x= Mxy(pos(find(fv==min(fv))),:)';
        pos=pos(find(fv==min(fv)));
        
    else //Equacoes de Igauladdae e Desigualdade 
     Aeqm=[Aeq(:,1:nv-m)  eye(ne(1)) -eye(ne(1))];
     Am=[A(:,1:nv-m) zeros(ni(1),2*ne(1))];
     lb=[lb(1:nv-m) zeros(1,2*ne(1))];
     ub=[ub(1:nv-m) Inf*ones(1,2*ne(1))];
     c=[f(1:nv-m) ones(1,2*ne(1))];
     op=options;
     fv=[];pos=[];Mxy=[];nlp=0;
     for i=1:2^ny
            Bm=B-f(nv-m+1:end)*M(i,:)';
            Beqm=Beq-f(nv-m+1:end)*M(i,:)';
            [X,FVAL,EXITFLAG,OUTPUT,LAMBDA]=linprog(c,Am,Bm,Aeqm,Beqm,lb,ub,x0,op);
            exitflag{i}=EXITFLAG;
            output{i}=OUTPUT;
            lambda{i}=LAMBDA;
            nlp=nlp+(EXITFLAG>0);
            xmy=[X;M(i,:)'];
            Mxy=[Mxy;xmy'];
            FV(i)=FVAL;
            FVALy(i)=c(1:nv-m)*X(1:nv-m)+f(nv-m+1:end)*M(i,:)';
            flagv(i)= sum((A*[xmy(1:nv-m);M(i,:)'] -B)>=0) + sum((Aeq*[xmy(1:nv-m);M(i,:)'] -Beq)>=tol);
            if EXITFLAG==1 
                if flagv(i)==0
                   fv=[fv;f*[xmy(1:nv-m);M(i,:)']];
                   pos=[pos;i];
                end
            end
        end
        fval=min(fv);
        MXY=[Mxy(:,1:nv-m) M];
        x=MXY(pos(find(fv==min(fv))),:)';
        pos=pos(find(fv==min(fv)));
    end
end
disp('====================================');
disp('====================================');
disp('RESULTADOS OBTIDOS:')
endfunction
