function [wo,St,x0,y0] = minlp(NLP,GRAD,n,m,rx,fig,x0,y0,xlb,xub,M,epsR,epsI,epsC,epsZ)
//     * * * Solution of MINLP * * *
//   Z = min(y,x) F(x,y) = C'y + f(x)
//    s.t. G(x,y) = By + g(x) <= 0
//         x in X in R^n
//         y in Y = {0,1}^m
//
//  usage: [wo,St,x0,y0] = minlp(NLP,GRAD,n,m,rx,fig,x0,y0,xlb,xub,M,epsR,epsI,epsC,epsZ)
//
//  wo      : best solution found [x,y,u,Z]
//            where u is the maximum constraint violation
//  St      : matrix = (number of constraints, number of LPs) for each NLP
//  NLP     : M file describing the NLP: [F,G] = nlp(x,y)
//            with the first line: if size(x,2) > n,  y = x(n+1:n+m); end
//  GRAD    : M file for gradients: [Fx,Gx,C,B] = grad(x,y)
//            with the first line: if size(x,2) > n,  y = x(n+1:n+m); end
//            and last lines: if size(x,2) > n, Fx = [Fx C]; Gx = [Gx B]; end
//                            Gx = Gx';
//  rx      : = 1 : relaxed initial NLP. = 0 : unrelaxed initial NLP (default)
//  fig     : = 1 : plot constraints surfaces. = 0 : no plots (default).
//  x0      : initial guess for x (default = zeros(1,n)). Returns first solution.
//  y0      : initial guess for y (default = round(rand(1,m)))
//  xlb     : lower bound for x (default = zeros(1,n))
//  xub     : upper bound for x (default = inf * ones(1,n))
//  M > 0   : constraint violation penalization (default = 1e3)
//  epsR    : real variable tolerance (default = 1e-4)
//  epsI    : integer variable tolerance (default = 1e-4)
//  epsC    : constraint violation tolerance (default = 1e-6)
//  epsZ    : objective function tolerance (default = 1e-4)
//

// authors: Argimiro R. Secchi and Elaine C. Pereira
//          Chemical Engineering Department
//          Applied Math Department
//          Universidade Federal do Rio Grande do Sul, Brasil
//          arge@enq.ufrgs.br
//
// version 1.0: 15/mar/1999

if nargin < 4
  disp('minimal usage: minlp(NLP,GRAD,n,m)');
  return;
end
if nargin < 5
  rx = 0;
end
if nargin < 6
  fig = 0;
end
if nargin < 7
  x0 = zeros(1,n);
end
if nargin < 8
  y0 = round(rand(1,m));
end
if nargin < 9
  xlb = zeros(1,n);
end
if nargin < 10
  xub = inf * ones(1,n);
end
if nargin < 11
  M = 1e3;
end
if nargin < 12
  epsR = 1e-4;
end
if nargin < 13
  epsI = 1e-4;
end
if nargin < 14
  epsC = 1e-6;
end
if nargin < 15
  epsZ = 1e-4;
end

if isempty(rx)
  rx = 0;
end
if isempty(fig)
  fig = 0;
end
if isempty(x0)
  x0 = zeros(1,n);
end
if isempty(y0)
  y0 = round(rand(1,m));
end
if isempty(xlb)
  xlb = zeros(1,n);
end
if isempty(xub)
  xub = inf * ones(1,n);
end
if isempty(M) | M <= 0
  M = 1e3;
end
if isempty(epsR)
  epsR = 1e-4;
end
if isempty(epsI)
  epsI = 1e-4;
end
if isempty(epsC)
  epsC = 1e-6;
end
if isempty(epsZ)
  epsZ = 1e-4;
end

ylb = zeros (1,m);
yub = ones(1,m);

//   * * * Solution of initial NLP * * *
//        Z = min(x) C'y0 + f(x) + Mu
//        s.t. By0 + g(x) <= u
//             x in X
//             y0 fixed
//             u in R+

if rx           // NLP with relaxed integrality
  x0=[x0 y0];
  wlb=[xlb ylb];
  wub=[xub yub];
else
  wlb=xlb;
  wub=xub;
end

op = foptions;
op(1) = 1;        // 1 - display some results; 0 - nothing
op(2) = epsR;
op(3) = epsZ;
op(4) = epsC;
op(13) = 0;
op(14) = 150;  // max function evaluation

// constr returns the best solution, even violating the constraints
// then u does not need be computed explicitly.

[x0,op,lambda] = constr(NLP,x0,op,wlb,wub,GRAD,y0);

if rx
  y0 = round(x0(n+1:n+m));
  dym = max(abs(x0(n+1:n+m)-y0));
  x0 = x0(1:n);
  if dym > epsI
    [x0,op,lambda] = constr(NLP,x0,op,xlb,xub,GRAD,y0);
  end   
end

// getting the values of f, g, and their gradients
[f_nlp,g_nlp] = evstr(NLP+"(x0,y0)");
u = max(g_nlp);             // u > 0 means constraint violation
[gf_nlp,gg_nlp,C,B] = evstr(GRAD+"(x0,y0)");
f = f_nlp - C * y0';
g = g_nlp' - B * y0';
gf = gf_nlp;
gg = gg_nlp';

if u > epsC                 // Z upper bound
  zu = inf;
else
  zu = f_nlp;
end

r = size(B,1);             // number of constraints
St = [r+1,0];              // matrix of NLPs: number of constraints x LPs

wo = [x0,y0,u,zu];         // best solution at each iteration

A = [gg,B,-ones(r,1),zeros(r,1); gf,C,M,-1];
b = [gg*x0'- g; gf*x0'-f]; 

wlb = [xlb,ylb,0,-inf];
wub = [xub,yub,inf,inf];

if fig, planos(n,A,b,y0,u,xlb,xub,1); end

Tr = [0, 1, 0, 0, 1:m];  // binary tree [parent, child, mj, nfix, fix, ~fix]
W = wo;                  // solution of the LPs
no = 1;                  // number of nodes in the tree

[flag,wo,A,b,W,Tr,no,St] = lp_nlp(NLP,GRAD,n,m,fig,wo,A,b,M,wlb,wub,W,Tr,1,no,epsR,epsI,epsC,epsZ,St);

while no
  if flag             // non integer solution
    w = W(no,:);      // create two child problems deleting the parent
    j = Tr(no,3);
    pai = Tr(no,2);
    nfx = Tr(no,4);
    fx = [Tr(no,5:4+nfx),j];
    vr = Tr(no,5+nfx:4+m);
    vr = vr(find(vr-j));
    Tr(no,:) = [pai,pai+1,0,nfx+1,fx,vr];
    no = no + 1;
    W(no,:) = w;
    Tr(no,:) = Tr(no-1,:);
    Tr(no,2) = pai+2;
    j = j + n;
    yj = round(w(j));
    W(no-1,j) = yj;
    W(no,j) = 1 - yj;
  end
  
// round trip for version 4.2
// version 5.2:   if no > 1 & Tr(no,1) == Tr(no-1,1)

  if no > 1
    xx = Tr(no-1,1);
  else
    xx = -1;
  end
  
  if Tr(no,1) == xx     // last two child problems with same parent
    old_no = no;
    [flag,wo,A,b,W,Tr,no,St] = lp_nlp(NLP,GRAD,n,m,fig,wo,A,b,M,wlb,wub,W,Tr,no,no,epsR,epsI,epsC,epsZ,St);
    if flag
      k = old_no - 1;
      [flag,wo,A,b,W,Tr,no,St] = lp_nlp(NLP,GRAD,n,m,fig,wo,A,b,M,wlb,wub,W,Tr,k,no,epsR,epsI,epsC,epsZ,St);
       if flag & old_no == no
         if W(old_no,n+m+2) > W(k,n+m+2)
           W([k old_no],:) = W([old_no k],:); // switch last two nodes
           Tr([k old_no],:) = Tr([old_no k],:);
         end
      end
    end
  else
    [flag,wo,A,b,W,Tr,no,St] = lp_nlp(NLP,GRAD,n,m,fig,wo,A,b,M,wlb,wub,W,Tr,no,no,epsR,epsI,epsC,epsZ,St);
  end
end

endfunction
