Programación estocástica

La formulación del problema de dos etapas anterior supone que los datos de la segunda etapa${\displaystyle \xi }$se modela como un vector aleatorio con una distribución de probabilidad conocida . Esto estaría justificado en muchas situaciones. Por ejemplo, la distribución de${\displaystyle \xi }$podría inferirse a partir de datos históricos si se supone que la distribución no cambia significativamente durante el período de tiempo considerado. Además, la distribución empírica de la muestra podría usarse como una aproximación a la distribución de los valores futuros de${\displaystyle \xi }$. Si uno tiene un modelo previo para${\displaystyle \xi }$, se podría obtener una distribución a posteriori mediante una actualización bayesiana.

Para resolver numéricamente el problema estocástico de dos etapas, a menudo es necesario suponer que el vector aleatorio${\displaystyle \xi }$tiene un número finito de posibles realizaciones, llamadas escenarios , por ejemplo${\displaystyle \xi _{1},\dots ,\xi _{K}}$, con sus respectivas masas de probabilidad${\displaystyle p_{1},\dots ,p_{K}}$Entonces, la esperanza en la función objetivo del problema de la primera etapa se puede escribir como la sumatoria: $${\displaystyle E[Q(x,\xi )]=\sum \limits _{k=1}^{K}p_{k}Q(x,\xi _{k})}$$ y, además, el problema de dos etapas puede formularse como un gran problema de programación lineal (esto se denomina el equivalente determinista del problema original, véase la sección §  Equivalente determinista de un problema estocástico ).

Cuando${\displaystyle \xi }$tiene un número infinito (o muy grande) de realizaciones posibles, el enfoque estándar es entonces representar esta distribución mediante escenarios. Este enfoque plantea tres preguntas, a saber:

1. Para saber cómo construir escenarios, consulte el apartado §  Construcción de escenarios ;
2. Cómo resolver el equivalente determinista. Los optimizadores como CPLEX y GLPK pueden resolver grandes problemas lineales/no lineales. El servidor NEOS, ^{[ 6 ]} alojado en la Universidad de Wisconsin, Madison , permite el acceso gratuito a muchos solucionadores modernos. La estructura de un equivalente determinista es particularmente adecuada para aplicar métodos de descomposición, ^{[ 7 ]} como la descomposición de Benders o la descomposición de escenarios;
3. Cómo medir la calidad de la solución obtenida con respecto al óptimo "verdadero".

Estas cuestiones no son independientes. Por ejemplo, el número de escenarios construidos afectará tanto a la viabilidad del equivalente determinista como a la calidad de las soluciones obtenidas.

Un programa lineal estocástico es una instancia específica del programa estocástico clásico de dos etapas. Un programa lineal estocástico se construye a partir de una colección de programas lineales (PL) multiperíodo, cada uno con la misma estructura pero con datos algo diferentes.${\displaystyle k^{th}}$LP de dos períodos, que representa el${\displaystyle k^{th}}$Este escenario puede considerarse que tiene la siguiente forma:

${\displaystyle {\begin{array}{lccccccc}{\text{Minimize}}&f^{T}x&+&g^{T}y&+&h_{k}^{T}z_{k}&&\\{\text{subject to}}&Tx&+&Uy&&&=&r\\&&&V_{k}y&+&W_{k}z_{k}&=&s_{k}\\&x&,&y&,&z_{k}&\geq &0\end{array}}}$

Los vectores${\displaystyle x}$y${\displaystyle y}$Contiene las variables del primer período, cuyos valores deben elegirse inmediatamente. El vector${\displaystyle z_{k}}$Contiene todas las variables para los períodos subsiguientes. Las restricciones${\displaystyle Tx+Uy=r}$Las restricciones incluyen únicamente variables del primer período y son las mismas en todos los escenarios. Las demás restricciones incluyen variables de períodos posteriores y difieren en algunos aspectos de un escenario a otro, lo que refleja la incertidumbre sobre el futuro.

Tenga en cuenta que resolver el${\displaystyle k^{th}}$La programación lineal de dos períodos es equivalente a suponer que${\displaystyle k^{th}}$escenario en el segundo período sin incertidumbre. Para incorporar incertidumbres en la segunda etapa, se deben asignar probabilidades a diferentes escenarios y resolver el equivalente determinista correspondiente.

En la práctica, podría ser posible construir escenarios recabando opiniones de expertos sobre el futuro. El número de escenarios construidos debería ser relativamente modesto para que el equivalente determinista obtenido pueda resolverse con un esfuerzo computacional razonable. A menudo se afirma que una solución óptima que utiliza solo unos pocos escenarios proporciona planes más adaptables que una que asume un solo escenario. En algunos casos, esta afirmación podría verificarse mediante una simulación. En teoría, existen algunas medidas que garantizan que una solución obtenida resuelva el problema original con una precisión razonable. Normalmente, en las aplicaciones, solo se considera la primera etapa de la solución óptima.${\displaystyle x^{*}}$tiene un valor práctico ya que casi siempre una realización "verdadera" de los datos aleatorios será diferente del conjunto de escenarios construidos (generados).

Suponer${\displaystyle \xi }$contiene${\displaystyle d}$componentes aleatorios independientes, cada uno de los cuales tiene tres posibles realizaciones (por ejemplo, las futuras realizaciones de cada parámetro aleatorio se clasifican como bajas, medias y altas), entonces el número total de escenarios es${\displaystyle K=3^{d}}$Tal crecimiento exponencial del número de escenarios hace que el desarrollo de modelos utilizando la opinión de expertos sea muy difícil incluso para un tamaño razonable.${\displaystyle d}$. La situación se vuelve aún peor si algunos componentes aleatorios de${\displaystyle \xi }$tienen distribuciones continuas.

Consideremos el siguiente problema de programación estocástica.

${\displaystyle \min \limits _{x\in X}\{g(x)=f(x)+E[Q(x,\xi )]\}}$

Aquí${\displaystyle X}$es un subconjunto cerrado no vacío de${\displaystyle \mathbb {R} ^{n}}$,${\displaystyle \xi }$es un vector aleatorio cuya distribución de probabilidad${\displaystyle P}$es compatible con un conjunto${\displaystyle \Xi \subset \mathbb {R} ^{d}}$, y${\displaystyle Q:X\times \Xi \rightarrow \mathbb {R} }$En el marco de la programación estocástica de dos etapas,${\displaystyle Q(x,\xi )}$viene dado por el valor óptimo del problema de segunda etapa correspondiente.

Supongamos que${\displaystyle g(x)}$está bien definido y tiene un valor finito para todos.${\displaystyle x\in X}$. Esto implica que para cada${\displaystyle x\in X}$el valor${\displaystyle Q(x,\xi )}$es finito casi con seguridad .

Supongamos que tenemos una muestra${\displaystyle \xi ^{1},\dots ,\xi ^{N}}$de${\displaystyle N}$realizaciones del vector aleatorio${\displaystyle \xi }$Esta muestra aleatoria puede considerarse como datos históricos de${\displaystyle N}$observaciones de${\displaystyle \xi }$o puede generarse mediante técnicas de muestreo de Monte Carlo. Entonces podemos formular una aproximación promedio de muestra correspondiente.

${\displaystyle \min \limits _{x\in X}\{{\hat {g}}_{N}(x)=f(x)+{\frac {1}{N}}\sum _{j=1}^{N}Q(x,\xi ^{j})\}}$

Por la ley de los grandes números tenemos que, bajo ciertas condiciones de regularidad${\displaystyle {\frac {1}{N}}\sum _{j=1}^{N}Q(x,\xi ^{j})}$converge puntualmente con probabilidad de 1 a${\displaystyle E[Q(x,\xi )]}$como${\displaystyle N\rightarrow \infty }$Además, bajo condiciones adicionales leves, la convergencia es uniforme. También tenemos${\displaystyle E[{\hat {g}}_{N}(x)]=g(x)}$, es decir,${\displaystyle {\hat {g}}_{N}(x)}$es un estimador insesgado de${\displaystyle g(x)}$Por lo tanto, es natural esperar que el valor óptimo y las soluciones óptimas del problema SAA converjan a sus contrapartes del problema verdadero como${\displaystyle N\rightarrow \infty }$.

Consistencia de los estimadores SAA

Supongamos que el conjunto factible${\displaystyle X}$del problema SAA es fijo, es decir, es independiente de la muestra. Sea${\displaystyle \vartheta ^{*}}$y${\displaystyle S^{*}}$sea ​​el valor óptimo y el conjunto de soluciones óptimas, respectivamente, del problema verdadero y sea${\displaystyle {\hat {\vartheta }}_{N}}$y${\displaystyle {\hat {S}}_{N}}$sean, respectivamente, el valor óptimo y el conjunto de soluciones óptimas del problema SAA.

1. Dejar${\displaystyle g:X\rightarrow \mathbb {R} }$y${\displaystyle {\hat {g}}_{N}:X\rightarrow \mathbb {R} }$Sea una secuencia de funciones (deterministas) de valor real. Las dos propiedades siguientes son equivalentes:
   • para cualquier${\displaystyle {\overline {x}}\in X}$y cualquier secuencia${\displaystyle \{x_{N}\}\subset X}$convergiendo a${\displaystyle {\overline {x}}}$resulta que${\displaystyle {\hat {g}}_{N}(x_{N})}$converge a${\displaystyle g({\overline {x}})}$
   • la función${\displaystyle g(\cdot )}$es continuo en${\displaystyle X}$y${\displaystyle {\hat {g}}_{N}(\cdot )}$converge a${\displaystyle g(\cdot )}$uniformemente en cualquier subconjunto compacto de${\displaystyle X}$
2. Si el objetivo del problema SAA${\displaystyle {\hat {g}}_{N}(x)}$converge al objetivo del problema verdadero${\displaystyle g(x)}$con probabilidad 1, como${\displaystyle N\rightarrow \infty }$, uniformemente en el conjunto factible${\displaystyle X}$. Entonces${\displaystyle {\hat {\vartheta }}_{N}}$converge a${\displaystyle \vartheta ^{*}}$con probabilidad 1 como${\displaystyle N\rightarrow \infty }$.
3. Supongamos que existe un conjunto compacto${\displaystyle C\subset \mathbb {R} ^{n}}$de tal manera que
   • el conjunto${\displaystyle S}$de soluciones óptimas del problema verdadero no es vacío y está contenido en${\displaystyle C}$
   • la función${\displaystyle g(x)}$es de valor finito y continuo en${\displaystyle C}$
   • la secuencia de funciones${\displaystyle {\hat {g}}_{N}(x)}$converge a${\displaystyle g(x)}$con probabilidad 1, como${\displaystyle N\rightarrow \infty }$, uniformemente en${\displaystyle x\in C}$
   • para${\displaystyle N}$lo suficientemente grande el conjunto${\displaystyle {\hat {S}}_{N}}$no está vacío y${\displaystyle {\hat {S}}_{N}\subset C}$con probabilidad 1

entonces${\displaystyle {\hat {\vartheta }}_{N}\rightarrow \vartheta ^{*}}$y${\displaystyle \mathbb {D} (S^{*},{\hat {S}}_{N})\rightarrow 0}$con probabilidad 1 como${\displaystyle N\rightarrow \infty }$. Tenga en cuenta que${\displaystyle \mathbb {D} (A,B)}$denota la desviación del conjunto${\displaystyle A}$del conjunto${\displaystyle B}$, definido como

${\displaystyle \mathbb {D} (A,B):=\sup _{x\in A}\{\inf _{x'\in B}\|x-x'\|\}}$

En algunas situaciones el conjunto factible${\displaystyle X}$del problema SAA se estima, luego el problema SAA correspondiente toma la forma

${\displaystyle \min _{x\in X_{N}}{\hat {g}}_{N}(x)}$

dónde${\displaystyle X_{N}}$es un subconjunto de${\displaystyle \mathbb {R} ^{n}}$dependiendo de la muestra y, por lo tanto, es aleatorio. Sin embargo, aún se pueden obtener resultados de consistencia para los estimadores SAA bajo algunos supuestos adicionales:

1. Supongamos que existe un conjunto compacto${\displaystyle C\subset \mathbb {R} ^{n}}$de tal manera que
   • el conjunto${\displaystyle S}$de soluciones óptimas del problema verdadero no es vacío y está contenido en${\displaystyle C}$
   • la función${\displaystyle g(x)}$es de valor finito y continuo en${\displaystyle C}$
   • la secuencia de funciones${\displaystyle {\hat {g}}_{N}(x)}$converge a${\displaystyle g(x)}$con probabilidad 1, como${\displaystyle N\rightarrow \infty }$, uniformemente en${\displaystyle x\in C}$
   • para${\displaystyle N}$lo suficientemente grande el conjunto${\displaystyle {\hat {S}}_{N}}$no está vacío y${\displaystyle {\hat {S}}_{N}\subset C}$con probabilidad 1
   • si${\displaystyle x_{N}\in X_{N}}$y${\displaystyle x_{N}}$converge con probabilidad 1 a un punto${\displaystyle x}$, entonces${\displaystyle x\in X}$
   • por algún momento${\displaystyle x\in S^{*}}$existe una secuencia${\displaystyle x_{N}\in X_{N}}$de tal manera que${\displaystyle x_{N}\rightarrow x}$con probabilidad 1.

entonces${\displaystyle {\hat {\vartheta }}_{N}\rightarrow \vartheta ^{*}}$y${\displaystyle \mathbb {D} (S^{*},{\hat {S}}_{N})\rightarrow 0}$con probabilidad 1 como${\displaystyle N\rightarrow \infty }$.

Asintótica del valor óptimo de SAA

Supongamos que la muestra${\displaystyle \xi ^{1},\dots ,\xi ^{N}}$es iid y fija un punto${\displaystyle x\in X}$. Luego, el estimador del promedio de la muestra${\displaystyle {\hat {g}}_{N}(x)}$, de${\displaystyle g(x)}$es imparcial y tiene varianza${\displaystyle {\frac {1}{N}}\sigma ^{2}(x)}$, dónde${\displaystyle \sigma ^{2}(x):=Var[Q(x,\xi )]}$Se supone que es finito. Además, por el teorema del límite central tenemos que

${\displaystyle {\sqrt {N}}[{\hat {g}}_{N}-g(x)]{\xrightarrow {\mathcal {D}}}Y_{x}}$

dónde${\displaystyle {\xrightarrow {\mathcal {D}}}}$denota convergencia en la distribución y${\displaystyle Y_{x}}$tiene una distribución normal con media${\displaystyle 0}$y varianza${\displaystyle \sigma ^{2}(x)}$, escrito como${\displaystyle {\mathcal {N}}(0,\sigma ^{2}(x))}$.

En otras palabras,${\displaystyle {\hat {g}}_{N}(x)}$tiene una distribución asintóticamente normal , es decir, para valores grandes${\displaystyle N}$,${\displaystyle {\hat {g}}_{N}(x)}$tiene aproximadamente una distribución normal con media${\displaystyle g(x)}$y varianza${\displaystyle {\frac {1}{N}}\sigma ^{2}(x)}$Esto conduce a lo siguiente (aproximado):${\displaystyle 100(1-\alpha )}$% intervalo de confianza para${\displaystyle f(x)}$:

${\displaystyle \left[{\hat {g}}_{N}(x)-z_{\alpha /2}{\frac {{\hat {\sigma }}(x)}{\sqrt {N}}},{\hat {g}}_{N}(x)+z_{\alpha /2}{\frac {{\hat {\sigma }}(x)}{\sqrt {N}}}\right]}$

dónde${\displaystyle z_{\alpha /2}:=\Phi ^{-1}(1-\alpha /2)}$(aquí${\displaystyle \Phi (\cdot )}$denota la función de distribución acumulada de la distribución normal estándar) y

${\displaystyle {\hat {\sigma }}^{2}(x):={\frac {1}{N-1}}\sum _{j=1}^{N}\left[Q(x,\xi ^{j})-{\frac {1}{N}}\sum _{j=1}^{N}Q(x,\xi ^{j})\right]^{2}}$

es la estimación de la varianza muestral de${\displaystyle \sigma ^{2}(x)}$. Es decir, el error de estimación de${\displaystyle g(x)}$es (estocásticamente) de orden${\displaystyle O({\sqrt {N}})}$.

La programación dinámica estocástica se utiliza frecuentemente para modelar el comportamiento animal en campos como la ecología del comportamiento . ^{[ 8 ] [ 9 ]} Las pruebas empíricas de modelos de forrajeo óptimo y transiciones del ciclo vital , como el abandono del nido en aves y la puesta de huevos en avispas parasitoides, han demostrado el valor de esta técnica de modelado para explicar la evolución de la toma de decisiones conductuales. Estos modelos suelen ser de varias etapas, en lugar de dos.

La programación dinámica estocástica es una herramienta útil para comprender la toma de decisiones en condiciones de incertidumbre. La acumulación de capital en condiciones de incertidumbre es un ejemplo; a menudo, los economistas de recursos la utilizan para analizar problemas bioeconómicos ^{[ 10 ]} donde interviene la incertidumbre, como el clima, etc.

El siguiente es un ejemplo de programación estocástica multietapa del ámbito financiero. Supongamos que en el momento${\displaystyle t=0}$Tenemos capital inicial${\displaystyle W_{0}}$invertir en${\displaystyle n}$activos. Supongamos además que se nos permite reequilibrar nuestra cartera en ocasiones.${\displaystyle t=1,\dots ,T-1}$pero sin inyectarle dinero adicional. En cada período${\displaystyle t}$Tomamos una decisión sobre la redistribución de la riqueza actual.${\displaystyle W_{t}}$entre los${\displaystyle n}$activos. Deje${\displaystyle x_{0}=(x_{10},\dots ,x_{n0})}$sean las cantidades iniciales invertidas en los n activos. Requerimos que cada${\displaystyle x_{i0}}$es no negativo y que la ecuación de balance${\displaystyle \sum _{i=1}^{n}x_{i0}=W_{0}}$debería mantenerse.

Considere los rendimientos totales${\displaystyle \xi _{t}=(\xi _{1t},\dots ,\xi _{nt})}$para cada período${\displaystyle t=1,\dots ,T}$Esto forma un proceso aleatorio con valores vectoriales.${\displaystyle \xi _{1},\dots ,\xi _{T}}$En ese período de tiempo${\displaystyle t=1}$, podemos reequilibrar la cartera especificando las cantidades${\displaystyle x_{1}=(x_{11},\dots ,x_{n1})}$invertidos en los respectivos activos. En ese momento se han realizado los rendimientos del primer período, por lo que es razonable utilizar esta información en la decisión de reequilibrio. Por lo tanto, las decisiones de la segunda etapa, en el momento${\displaystyle t=1}$, son en realidad funciones de realización del vector aleatorio${\displaystyle \xi _{1}}$, es decir,${\displaystyle x_{1}=x_{1}(\xi _{1})}$. De manera similar, en el momento${\displaystyle t}$la decisión${\displaystyle x_{t}=(x_{1t},\dots ,x_{nt})}$es una función${\displaystyle x_{t}=x_{t}(\xi _{[t]})}$de la información disponible proporcionada por${\displaystyle \xi _{[t]}=(\xi _{1},\dots ,\xi _{t})}$la historia del proceso aleatorio hasta el momento${\displaystyle t}$Una secuencia de funciones${\displaystyle x_{t}=x_{t}(\xi _{[t]})}$,${\displaystyle t=0,\dots ,T-1}$, con${\displaystyle x_{0}}$Al ser constante, define una política implementable del proceso de decisión. Se dice que dicha política es factible si satisface las restricciones del modelo con probabilidad 1, es decir, las restricciones de no negatividad.${\displaystyle x_{it}(\xi _{[t]})\geq 0}$,${\displaystyle i=1,\dots ,n}$,${\displaystyle t=0,\dots ,T-1}$y el equilibrio de las limitaciones de riqueza,

${\displaystyle \sum _{i=1}^{n}x_{it}(\xi _{[t]})=W_{t},}$

donde en el período${\displaystyle t=1,\dots ,T}$la riqueza${\displaystyle W_{t}}$es dado por

${\displaystyle W_{t}=\sum _{i=1}^{n}\xi _{it}x_{i,t-1}(\xi _{[t-1]}),}$

lo cual depende de la realización del proceso aleatorio y de las decisiones tomadas hasta el momento.${\displaystyle t}$.

Supongamos que el objetivo es maximizar la utilidad esperada de esta riqueza en el último período, es decir, considerar el problema

${\displaystyle \max E[U(W_{T})].}$

Este es un problema de programación estocástica multietapa, donde las etapas están numeradas desde${\displaystyle t=0}$a${\displaystyle t=T-1}$La optimización se realiza sobre todas las políticas implementables y factibles. Para completar la descripción del problema, también es necesario definir la distribución de probabilidad del proceso aleatorio.${\displaystyle \xi _{1},\dots ,\xi _{T}}$Esto se puede hacer de varias maneras. Por ejemplo, se puede construir un árbol de escenarios particular que defina la evolución temporal del proceso. Si en cada etapa se permite que el rendimiento aleatorio de cada activo tenga dos continuaciones, independientemente de otros activos, entonces el número total de escenarios es${\displaystyle 2^{nT}.}$

Para escribir ecuaciones de programación dinámica , considere el problema multietapa anterior hacia atrás en el tiempo. En la última etapa${\displaystyle t=T-1}$, una constatación${\displaystyle \xi _{[T-1]}=(\xi _{1},\dots ,\xi _{T-1})}$del proceso aleatorio es conocido y${\displaystyle x_{T-2}}$ha sido elegido. Por lo tanto, es necesario resolver el siguiente problema.

${\displaystyle {\begin{array}{lrclr}\max \limits _{x_{T-1}}&E[U(W_{T})|\xi _{[T-1]}]&\\{\text{subject to}}&W_{T}&=&\sum _{i=1}^{n}\xi _{iT}x_{i,T-1}\\&\sum _{i=1}^{n}x_{i,T-1}&=&W_{T-1}\\&x_{T-1}&\geq &0\end{array}}}$

dónde${\displaystyle E[U(W_{T})|\xi _{[T-1]}]}$denota la esperanza condicional de${\displaystyle U(W_{T})}$dado${\displaystyle \xi _{[T-1]}}$. El valor óptimo del problema anterior depende de${\displaystyle W_{T-1}}$y${\displaystyle \xi _{[T-1]}}$y se denota${\displaystyle Q_{T-1}(W_{T-1},\xi _{[T-1]})}$.

De manera similar, en las etapas${\displaystyle t=T-2,\dots ,1}$uno debería resolver el problema

${\displaystyle {\begin{array}{lrclr}\max \limits _{x_{t}}&E[Q_{t+1}(W_{t+1},\xi _{[t+1]})|\xi _{[t]}]&\\{\text{subject to}}&W_{t+1}&=&\sum _{i=1}^{n}\xi _{i,t+1}x_{i,t}\\&\sum _{i=1}^{n}x_{i,t}&=&W_{t}\\&x_{t}&\geq &0\end{array}}}$

cuyo valor óptimo se denota por${\displaystyle Q_{t}(W_{t},\xi _{[t]})}$. Finalmente, en la etapa${\displaystyle t=0}$uno resuelve el problema

${\displaystyle {\begin{array}{lrclr}\max \limits _{x_{0}}&E[Q_{1}(W_{1},\xi _{[1]})]&\\{\text{subject to}}&W_{1}&=&\sum _{i=1}^{n}\xi _{i,1}x_{i0}\\&\sum _{i=1}^{n}x_{i0}&=&W_{0}\\&x_{0}&\geq &0\end{array}}}$

Todos los problemas de programación estocástica discreta pueden representarse con cualquier lenguaje de modelado algebraico , implementando manualmente la no anticipación explícita o implícita para asegurar que el modelo resultante respete la estructura de la información disponible en cada etapa. Una instancia de un problema de programación estocástica generada por un lenguaje de modelado general tiende a crecer considerablemente (linealmente con el número de escenarios), y su matriz pierde la estructura intrínseca a esta clase de problemas, que de otro modo podría aprovecharse en el momento de la solución mediante algoritmos de descomposición específicos. Están empezando a aparecer extensiones a lenguajes de modelado diseñados específicamente para programación estocástica; véase:

• AIMMS – admite la definición de problemas SP
• EMP SP (Programación Matemática Extendida para Programación Estocástica): un módulo de GAMS creado para facilitar la programación estocástica (incluye palabras clave para distribuciones paramétricas, restricciones de probabilidad y medidas de riesgo como Valor en Riesgo y Déficit Esperado ).
• SAMPL : un conjunto de extensiones de AMPL diseñadas específicamente para expresar programas estocásticos (incluye sintaxis para restricciones de probabilidad, restricciones de probabilidad integradas y problemas de optimización robusta ).

Ambos pueden generar un formato a nivel de instancia SMPS, que transmite de forma no redundante la estructura del problema al solucionador.