function [k, X, Y, pasos, mu] = vnpls(F, signoX, signoXF, k, X)
%vnpls
%   [k, X, Y] = vnpls(F)
% 
%   Busca la solución de VN mediante la programación lineal secuencial.
%   Resuelve sucesivamente la aproximación lineal de VN
%      max k                                       min -mu
%      X F(kn) + k Xn F'(kn) - kn Xn F'(kn) ~ 0    -Y ~ 0
%      sum(X) - 1 = 0                              -mu >=< 0
%      X ~ 0                                       F(kn) Y + mu ~ 0
%      k >= 0                                      1 + Xn F'(kn) Y <= 0
%   hasta que las variables converjan.
% 
%   [k, X, Y] = vnpls(F, signoX, signoXF, k0, X0)
%      F         recetas en tiempo discreto o continuo
%      signoX    signo de las intensidades; por defecto X >= 0
%      signoXF   signo de los balances materiales; por defecto X F(k) = 0
%      k0        factor con el que se inicia la iteración
%      X0        intensidades con las que se inicia la iteración
%      k, X, Y   factor, intensidades y valores solución
% 
%   Los códigos para los signos son   0 =    1 >=    2 <=    3 >=<

%   necesita forma0.m, forma1.m, programalineal.m

if nargin < 2, signoX = []; end
if nargin < 3, signoXF = []; end
if nargin < 4, k = []; end
if nargin < 5, X = []; end

cerosalida = 1e-14;                                             %define la precisión del cero
if isempty(k), k = 1; end                                       %define kn como 1 si no lo está ya
Fn = forma0(F, k);                                              %calcula F(kn)
[filas, columnas] = size(Fn);                                   %determina el tamaño del problema
if isempty(X), X = ones(1, filas) ./ filas; end                 %define Xn si no lo está ya
C = [zeros(filas, 1); 1];                                       %define la función objetivo como 0 X + 1 k
Aeq = [ones(filas, 1); 0];                                      %define la normalización sum(X) - 1 = 0
Deq = -1;                                                       %
if ~isempty(signoX), Aeq(signoX==2) = -1; end                   %si X <= 0 adecua la normalización
for pasos = 1:50                                                %inicia una iteración de hasta 50 pasos
    memvar = [k, X];                                                %almacena las variables
    XDFn = X * forma1(F, k);                                        %calcula Xn F'(kn)
    A = [Fn; XDFn];                                                 %define |X, k| |F(kn)    | + |-kn Xn F'(kn)|
    D = -k * XDFn;                                                  %              |Xn F'(kn)|
    [U,Z]=programalineal([A,Aeq],C,[D,Deq],[signoX,1],[signoXF,0]); %resuelve el programa lineal
    X=U(1:filas); k=U(filas+1); Y=Z(1:columnas); mu=Z(columnas+1);  %extrae las variables de la solución
    if abs([k, X] - memvar) < cerosalida, return; end               %si las variables cambian poco interrumpe
    Fn = forma0(F, k);                                              %calcula F(kn)
end                                                             %repite la iteración
disp('Atención, se ha superado el máximo de iteraciones');      %avisa de que se ha superado el máximo de iteraciones

