function [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)

// 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

//   * * * Solution of the k-th LP * * *
//   zl = min(x,y,u,alpha)   alpha
//   s.t.   alfa >= C'y + f(xj) + gf(xj)'(x-xj) + Mu
//          By + g(xj) + gg(xj)'(x-xj) <= u
//          x in X               
//          y in Y
//          alpha in R+

nm = n+m;
zu = wo(nm+2);
w = W(k,:);
w(nm+1) = 0;        // u = 0 (initial guess for LP)
if w(nm+2) == inf
  w(nm+2) = 0;      // alpha = 0 (initial guess for LP)
end

nfx = Tr(k,4);
nlps = size(St,1);

if nfx
  fx = n+Tr(k,5:4+nfx);                     // yj fixed in w
  vr = [1:n,n+Tr(k,5+nfx:4+m),nm+1,nm+2];   // w - yj fixed
  br = b - A(:,fx)*w(fx)';
  h = [zeros(nm-nfx+1,1); 1];
  wr = linprog(h,A(:,vr),br,[],[],wlb(vr),wub(vr),w(vr))';
  //wr = lp(h,A(:,vr),br,wlb(vr),wub(vr),w(vr))';
  w(vr) = wr;
else
  h = [zeros(nm+1,1); 1];
  w = linprog(h,A,b,[],[],wlb,wub,w)';
  //w = lp(h,A,b,wlb,wub,w)';
end

St(nlps,2) = St(nlps,2) + 1;              // LPs counter

//if w(nm+1) > epsC     // constraint violation
//  w(nm+2) = inf;
//end

W(k,:) = w;
zl = w(nm+2);         // lower bound for Z

if zl+epsZ >= zu           // deleting k-th node of the tree
  flag = 1;
  W = W([1:k-1, k+1:no],:);
  Tr = Tr([1:k-1, k+1:no],:);  // delete all P(i,k) such that z(i,k) >= zu
  no = no - 1;
  return;
end

x0 = w(1:n);
y = w(n+1:n+m);
[dym,j] = max(abs(y-round(y)));

if dym < epsI
  flag = 0;        // integer solution of last LP
  y = round(y);
  Tr(no,3) = 0;

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

  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.

  [x,op,lambda] = constr(NLP,x0,op,wlb(1:n),wub(1:n),GRAD,y);

// getting the values of f, g, and their gradients
  [f_nlp,g_nlp] = eval([NLP '(x,y)']);
  u = max(g_nlp);
  [gf_nlp,gg_nlp,C,B] = eval([GRAD '(x,y)']);
  f = f_nlp - C * y';
  g = g_nlp' - B * y';
  gf = gf_nlp;
  gg = gg_nlp';
  
  if u > epsC
    znl = inf;
    lambda=-lambda(1:size(B,1));
//   A = [A; zeros(1,n),C+lambda'*B,M,-1];   // adding Bender's cut
//   b = [b; -f-lambda'*g];
//   A = [A; zeros(1,n),lambda'*B,-sum(lambda),0];   // adding Bender's cut
//   b = [b; -lambda'*g];
//   A = [A; zeros(1,n),lambda'*B,0,0];   // adding Bender's cut
//   b = [b; -lambda'*g];
    A = [A; zeros(1,n),C+lambda'*B,M-sum(lambda),-1];   // adding Bender's cut
    b = [b; -f-lambda'*g];
  else
    znl = f_nlp;
    A = [A; gf,C,M,-1];                     // adding OA
    b = [b; gf*x'-f];
  end
  
  if fig, planos(n,A,b,y,u,wlb(1:n),wub(1:n),1); end
  
  if znl < zu
    zu = znl;
    wo = [x,y,u,zu];
  end
  
  nlps = nlps + 1;                // NLPs counter
  St(nlps,1:2) = [size(A,1), 0];

// round trip for version 4.2
// version 5.2:   while no > 0 & W(no,nm+2) >= zu

  if no > 0
    xx = W(no,nm+2);
  else
    xx = -inf;
  end
  while no > 0 && xx >= zu
       no = no - 1;
       W = W(1:no,:);
       Tr = Tr(1:no,:);  // delete all P(i,k) such that z(i,k) >= zu
       if no > 0, xx = W(no,nm+2); end
  end
else
  Tr(k,3) = j;    // less integer yj
  flag = 1;       // non-integer solution
end

endfunction

