vnpls.m

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