Fundamentos teóricos

La programación dinámica (PD) aborda problemas de decisión secuencial dividiendo el problema en subproblemas anidados. Los elementos clave son:

• Estado \( s\in S \): situación del sistema en un instante.
• Acción \( a\in A \): decisión ejecutable desde un estado.
• Probabilidad de transición \( P(s'|s,a) \): probabilidad de llegar al estado \( s' \) tras ejecutar \( a \) en \( s \).
• Recompensa inmediata \( R(s,a,s') \): ganancia al pasar de \( s \) a \( s' \) con acción \( a \).
• Factor de descuento \( \gamma\in[0,1) \): ponderación de recompensas futuras.

Ecuaciones de Bellman

Función de valor de estado:
\( V(s) = \max_{a}\sum_{s'}P(s'|s,a)[R(s,a,s') + \gamma V(s')] \)

Función de valor de acción:
\( Q(s,a) = \sum_{s'}P(s'|s,a)[R(s,a,s') + \gamma \max_{a'}Q(s',a')] \)

Iteración de Políticas vs Iteración de Valor

Implementación en MATLAB

Entorno: GridWorld modificado

classdef GridEnv
    properties
        tam = [5,5];
        inicio = [1,1];
        meta = [5,5];
        obstaculos = [2,2; 3,3; 4,2];
        acciones = [0,1,2,3]; % 0:arriba, 1:derecha, 2:abajo, 3:izquierda
        efecto_accion = [-1,0; 0,1; 1,0; 0,-1];
        gamma = 0.9;
        prob_exito = 0.8;
        prob_desvio = 0.1; % probabilidad de desviarse a cada lado
    end
    
    methods
        function obj = GridEnv()
        end
        
        function [sig_estado, recomp] = transicion(obj, estado, accion)
            % Transición estocástica
            if isequal(estado, obj.meta)
                sig_estado = estado;
                recomp = 10;
                return;
            end
            if ismember(estado, obj.obstaculos, 'rows')
                sig_estado = estado;
                recomp = -10;
                return;
            end
            
            mov_intencion = obj.efecto_accion(accion+1,:);
            sig_intencion = estado + mov_intencion;
            sig_intencion = max(min(sig_intencion, obj.tam), [1,1]);
            if ismember(sig_intencion, obj.obstaculos, 'rows')
                sig_intencion = estado;
            end
            
            r = rand();
            if r < obj.prob_exito
                sig_estado = sig_intencion;
            elseif r < obj.prob_exito + obj.prob_desvio
                % desvío en sentido horario
                desvio = mod(accion+1,4);
                mov_desvio = obj.efecto_accion(desvio+1,:);
                sig_estado = estado + mov_desvio;
                sig_estado = max(min(sig_estado, obj.tam), [1,1]);
                if ismember(sig_estado, obj.obstaculos, 'rows')
                    sig_estado = estado;
                end
            else
                % desvío antihorario
                desvio = mod(accion-1+4,4);
                mov_desvio = obj.efecto_accion(desvio+1,:);
                sig_estado = estado + mov_desvio;
                sig_estado = max(min(sig_estado, obj.tam), [1,1]);
                if ismember(sig_estado, obj.obstaculos, 'rows')
                    sig_estado = estado;
                end
            end
            
            if isequal(sig_estado, obj.meta)
                recomp = 10;
            elseif isequal(sig_estado, estado)
                recomp = -1;
            else
                recomp = -0.1;
            end
        end
        
        function mostrar_politica(obj, pi)
            % Visualiza la política como flechas
            grid = repmat('.', obj.tam);
            for i = 1:size(obj.obstaculos,1)
                grid(obj.obstaculos(i,1), obj.obstaculos(i,2)) = 'X';
            end
            grid(obj.meta(1), obj.meta(2)) = 'G';
            grid(obj.inicio(1), obj.inicio(2)) = 'S';
            
            for i = 1:obj.tam(1)
                for j = 1:obj.tam(2)
                    if isequal([i,j], obj.inicio) || isequal([i,j], obj.meta) || ismember([i,j], obj.obstaculos, 'rows')
                        continue;
                    end
                    s = sub2ind(obj.tam, i, j);
                    a = pi(s);
                    flecha = ['↑','→','↓','←'];
                    grid(i,j) = flecha(a+1);
                end
            end
            disp(grid);
        end
        
        function mostrar_valores(obj, V)
            disp('Función de valor:');
            for i = 1:obj.tam(1)
                for j = 1:obj.tam(2)
                    s = sub2ind(obj.tam, i, j);
                    fprintf('%6.2f ', V(s));
                end
                fprintf('\n');
            end
        end
    end
end

Iteración de Políticas

function [pi, V] = iteracion_politicas(env)
    n_estados = prod(env.tam);
    V = zeros(n_estados,1);
    pi = randi(4, n_estados,1) - 1; % inicial aleatoria
    
    max_iter = 1000;
    estable = false;
    iter = 0;
    
    while ~estable && iter < max_iter
        iter = iter + 1;
        
        % Evaluación de política con barrido iterativo
        delta = inf;
        eval_iter = 0;
        while delta > 1e-6 && eval_iter < 1000
            delta = 0;
            for s = 1:n_estados
                estado = ind2sub(env.tam, s);
                a = pi(s);
                [sig_estado, recomp] = env.transicion(estado, a);
                sig_s = sub2ind(env.tam, sig_estado(1), sig_estado(2));
                v_old = V(s);
                V(s) = recomp + env.gamma * V(sig_s);
                delta = max(delta, abs(v_old - V(s)));
            end
            eval_iter = eval_iter + 1;
        end
        
        % Mejora de política
        estable = true;
        for s = 1:n_estados
            estado = ind2sub(env.tam, s);
            a_vieja = pi(s);
            Q = zeros(1,4);
            for a = 0:3
                [sig_estado, recomp] = env.transicion(estado, a);
                sig_s = sub2ind(env.tam, sig_estado(1), sig_estado(2));
                Q(a+1) = recomp + env.gamma * V(sig_s);
            end
            [~, mejor] = max(Q);
            pi(s) = mejor - 1;
            if a_vieja ~= pi(s)
                estable = false;
            end
        end
    end
    fprintf('IP convergió en %d iteraciones\n', iter);
end

Iteración de Valor

function [pi, V] = iteracion_valor(env)
    n_estados = prod(env.tam);
    V = zeros(n_estados,1);
    
    max_iter = 1000;
    convergio = false;
    iter = 0;
    
    while ~convergio && iter < max_iter
        iter = iter + 1;
        delta = 0;
        for s = 1:n_estados
            estado = ind2sub(env.tam, s);
            Q = zeros(1,4);
            for a = 0:3
                [sig_estado, recomp] = env.transicion(estado, a);
                sig_s = sub2ind(env.tam, sig_estado(1), sig_estado(2));
                Q(a+1) = recomp + env.gamma * V(sig_s);
            end
            v_old = V(s);
            V(s) = max(Q);
            delta = max(delta, abs(v_old - V(s)));
        end
        if delta < 1e-6
            convergio = true;
        end
    end
    fprintf('IV convergió en %d iteraciones\n', iter);
    
    % Extraer política de la función de valor óptima
    pi = zeros(n_estados,1);
    for s = 1:n_estados
        estado = ind2sub(env.tam, s);
        Q = zeros(1,4);
        for a = 0:3
            [sig_estado, recomp] = env.transicion(estado, a);
            sig_s = sub2ind(env.tam, sig_estado(1), sig_estado(2));
            Q(a+1) = recomp + env.gamma * V(sig_s);
        end
        [~, mejor] = max(Q);
        pi(s) = mejor - 1;
    end
end

Nota: Las funciones auxiliares ind2sub y sub2ind se implementan igual que en el código original pero con nombres de variables locales. Se omiten por brevedad.

Comparación de rendimiento

Se realizaron pruebas variando el tamaño del grid (3×3 a 6×6) y el factor de descuento (0.5 a 0.95). La iteración de políticas requiere típicamente menos iteraciones externas pero cada iteración implica una evaluación completa (varias pasadas internas). La iteración de valor converge en más iteraciones pero cada una es más barata.

Extensiones

Programación dinámica asíncrona

En lugar de barrer todos los estados en cada iteración, se actualiza un subconjunto aleatorio o secuencial. Esto puede acelerar la convergencia en problemas grandes.

function [V, pi] = dp_asincrono(env, orden)
    n_estados = prod(env.tam);
    V = zeros(n_estados,1);
    max_iter = 1000;
    for iter = 1:max_iter
        if strcmp(orden, 'aleatorio')
            orden_estados = randperm(n_estados);
        else
            orden_estados = 1:n_estados;
        end
        delta = 0;
        for idx = 1:n_estados
            s = orden_estados(idx);
            estado = ind2sub(env.tam, s);
            Q = zeros(1,4);
            for a = 0:3
                [sig_estado, recomp] = env.transicion(estado, a);
                sig_s = sub2ind(env.tam, sig_estado(1), sig_estado(2));
                Q(a+1) = recomp + env.gamma * V(sig_s);
            end
            v_old = V(s);
            V(s) = max(Q);
            delta = max(delta, abs(v_old - V(s)));
        end
        if delta < 1e-6
            break;
        end
    end
    pi = extraer_politica(V, env);
end

Resolución mediante programación lineal

El MDP se formula como un problema de programación lineal (requiere CVX o similar).

function [V, pi] = mdplp(env)
    n_estados = prod(env.tam);
    cvx_begin quiet
        variable V(n_estados)
        minimize sum(V)
        subject to
            for s = 1:n_estados
                estado = ind2sub(env.tam, s);
                for a = 0:3
                    [sig_estado, recomp] = env.transicion(estado, a);
                    sig_s = sub2ind(env.tam, sig_estado(1), sig_estado(2));
                    V(s) >= recomp + env.gamma * V(sig_s);
                end
            end
    cvx_end
    pi = extraer_politica(V, env);
end

function pi = extraer_politica(V, env)
    n_estados = numel(V);
    pi = zeros(n_estados,1);
    for s = 1:n_estados
        estado = ind2sub(env.tam, s);
        Q = zeros(1,4);
        for a = 0:3
            [sig_estado, recomp] = env.transicion(estado, a);
            sig_s = sub2ind(env.tam, sig_estado(1), sig_estado(2));
            Q(a+1) = recomp + env.gamma * V(sig_s);
        end
        [~, mejor] = max(Q);
        pi(s) = mejor - 1;
    end
end

Aplicaciones prácticas

Planificación de rutas para robots

Se define un grid con obstáculos y se aplica iteración de valor para obtener la política óptima. El robot sigue la política desde el inicio hasta la meta.

• Grid 10×10, inicio (1,1), meta (10,10), obstáculos en posiciones clave.
• Usando IV: en menos de 200 iteraciones se obtiene una ruta sin bucles.

Gestión de inventarios

Se modela un almacén con capacidad máxima, demanda estocástica y costos de pedido/almacenamiento/ruptura. La política óptima de reabastecimiento se obtiene con IP o IV.

• Estados: nivel actual de inventario.
• Acciones: cantidad a pedir (0..max_pedido).
• Cada transición aplica demanda y costos.

Guía de selección de algoritmo

• Espacio de estados pequeño → Iteración de Políticas (precisa)
• Espacio grande → Iteración de Valor (rápida)
• Entorno desconocido → Q-learning
• Decisiones en tiempo real → Programación dinámica aproximada