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];
  if version > 4,
    wr = linprog(h,A(:,vr),br,[],[],wlb(vr),wub(vr),w(vr))';
  else
    wr = lp(h,A(:,vr),br,wlb(vr),wub(vr),w(vr))';
  end
  w(vr) = wr;
else
  h = [zeros(nm+1,1); 1];
  if version > 4,
    w = linprog(h,A,b,[],[],wlb,wub,w)';
  else
    w = lp(h,A,b,wlb,wub,w)';
  end
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

