Implementación del Método de Conjunto Activo y Simplex con Factorización QR en MATLAB
Implementación del Método de Conjunto Activo (Simplex QR)
El siguiente código presenta la implementación en MATLAB de una variante del método de conjunto activo (Active-Set Method) para problemas de programación lineal con cotas y restricciones generales, empleando actualizaciones de la factorización QR.
Parámetros de la Función
Entradas
- x: Vector inicial de variables de decisión.
- c: Vector de coeficientes de la función objetivo lineal.
- Aeq, beq: Matriz y vector de restricciones lineales de igualdad ($A_{eq} x = b_{eq}$).
- A, b: Matriz y vector de restricciones lineales de desigualdad ($A x \le b$).
- ub: Vector de cotas superiores ($x \le ub$).
- lb: Vector de cotas inferiores ($x \ge lb$).
- itermax: Número máximo de iteraciones permitidas.
Salidas
- x: Solución óptima calculada.
- fx: Valor de la función objetivo en el óptimo ($c^T x$).
- lambda: Multiplicadores de Lagrange asociados a las restricciones activas.
- itrab: Índices de las restricciones que conforman el conjunto de trabajo activo.
- iredes: Vector indicador del estado de las restricciones de desigualdad.
- iter: Número de iteraciones ejecutadas.
- ifallo: Indicador de convergencia o error (0 si convergió con éxito, 1 si falló).
Código Fuente Corregido
function [x, fx, lambda, itrab, iredes, iter, ifallo] = simplex7(x, c, Aeq, A, beq, b, ub, lb, itermax)
eps1 = eps^0.75;
fx = []; lambda = []; itrab = []; iredes = []; iter = 0; ifallo = 0;
% Calculamos el conjunto de trabajo inicial con las restricciones de igualdad
Ak = Aeq;
[Q, R, ifallo2] = qr(Ak');
if ifallo2 == 1
return;
end
% Determinación de restricciones de desigualdad activas
nd = length(A(:, 1)) + length(ub) + length(lb); % Número total de restricciones de desigualdad
iredes = zeros(nd, 1);
mi = length(Aeq(:, 1));
m = length(Aeq(:, 1));
n = length(x);
I = eye(n);
% Verificación de restricciones de desigualdad A*x <= b
for i = 1:length(A(:, 1))
if (A(i, :) * x - b(i)) >= -eps1
[Q, R, ind] = anfm01(Q, R, A(i, :)', n, m, 0);
if ind == 0
m = m + 1;
iredes(i) = 1;
itrab(m) = mi + i;
end
end
if m == n
break;
end
end
% Verificación de cotas inferiores lb <= x
for i = 1:n
if (lb(i) - x(i)) >= -eps1
[Q, R, ind] = anfm01(Q, R, -I(i, :)', n, m, -i);
if ind == 0
m = m + 1;
iredes(i + length(A(:, 1))) = 1;
itrab(m) = length(A(:, 1)) + mi + i;
end
end
if m == n
break;
end
end
% Verificación de cotas superiores x <= ub
for i = 1:n
if (ub(i) - x(i)) <= eps1
[Q, R, ind] = anfm01(Q, R, I(i, :)', n, m, i);
if ind == 0
m = m + 1;
iredes(i + length(A(:, 1)) + length(lb)) = 1;
itrab(m) = length(A(:, 1)) + length(lb) + mi + i;
end
end
if m == n
break;
end
end
% Bucle principal de iteraciones del Método de Conjunto Activo
for j = 1:itermax
Z = Q(:, m + 1:n);
if norm(Z' * c) <= eps1
[lambda, ifallo] = solvetrsup(R(1:m, 1:m), -Q(:, 1:m)' * c);
if ifallo == 1
fprintf('Se ha producido un error al resolver el sistema triangular.\n');
return;
end
[valor, icol] = min(lambda(mi + 1:m)); % Posición del multiplicador más negativo
if valor >= -eps1
fprintf('Todos los multiplicadores son no negativos. Se ha encontrado una solución: %13.5f\n', c' * x);
fx = c' * x;
return;
else
% Eliminación de restricción del conjunto de trabajo
[Q, R] = anfm02(Q, R, n, m, mi + icol);
m = m - 1;
iredes(itrab(mi + icol) - mi) = 0;
itrab(mi + icol) = [];
Z = Q(:, m + 1:n); % Matriz Z actualizada
d = (-sign(A(:, icol)' * Z(:, 1))) * Z(:, 1);
end
else
d = -Z * (Z' * c); % Dirección en el espacio nulo
end
% Cálculo del tamaño de paso ro (step length)
coord1 = find(iredes(1:length(A(:, 1))) == 0 & (A * d) > eps1);
ro = 1;
posicionro1 = 0; posicionro2 = 0; posicionro3 = 0;
for i = 1:length(coord1)
ro1 = (b(coord1(i)) - A(coord1(i), :) * x) / (A(coord1(i), :) * d);
if ro1 < ro
ro = ro1;
posicionro1 = i;
end
end
coord2 = find(iredes(length(A(:, 1)) + 1:length(A(:, 1)) + n) == 0 & -d > eps1);
for i = 1:length(coord2)
ro2 = (lb(coord2(i)) - x(coord2(i))) / (-d(coord2(i)));
if ro2 < ro
ro = ro2;
posicionro2 = i;
posicionro1 = 0;
end
end
coord3 = find(iredes(length(A(:, 1)) + n + 1:nd) == 0 & d > eps1);
for i = 1:length(coord3)
ro3 = (ub(coord3(i)) - x(coord3(i))) / (d(coord3(i)));
if ro3 < ro
ro = ro3;
posicionro3 = i;
posicionro1 = 0;
posicionro2 = 0;
end
end
% Actualización de la solución y adición de la nueva restricción activa
x = x + ro * d;
if posicionro1 ~= 0
[Q, R, ind] = anfm01(Q, R, A(coord1(posicionro1), :)', n, m, 0);
m = m + 1;
itrab(m) = mi + coord1(posicionro1);
iredes(coord1(posicionro1)) = 1;
elseif posicionro2 ~= 0
[Q, R, ind] = anfm01(Q, R, -I(coord2(posicionro2), :)', n, m, -coord2(posicionro2));
m = m + 1;
itrab(m) = mi + length(A(:, 1)) + coord2(posicionro2);
iredes(length(A(:, 1)) + coord2(posicionro2)) = 1;
elseif posicionro3 ~= 0
[Q, R, ind] = anfm01(Q, R, I(coord3(posicionro3), :)', n, m, coord3(posicionro3));
m = m + 1;
itrab(m) = mi + n + length(A(:, 1)) + coord3(posicionro3);
iredes(length(A(:, 1)) + n + coord3(posicionro3)) = 1;
end
iter = iter + 1;
end
endAspectos Matemáticos Relevantes
- Factorización QR y Espacio Nulo: Las columnas de $Q$ particionadas en $Z = Q(:, m+1:n)$ forman una base ortonormal del espacio nulo del conjunto de restricciones activas.
- Condiciones de Optimalidad KKT: Cuando la proyección en el gradiente reducido $Z^T c$ es prácticamente nula, se evalúan los multiplicadores de Lagrange mediante el sistema triangular $R^T \lambda = -Q_1^T c$.
- Criterio de Parada: Si todos los multiplicadores asociados a desigualdades satisfacen $\lambda_i \ge 0$, se garantiza la condición de mínimo local/global.