Planejamento da produção sob incerteza: programação estocástica versus otimização robusta

Production planning under uncertainty: stochastic programming versus robust optimization

Resumos

Otimizar problemas de planejamento da produção sob incertezas é um desafio, pois é preciso definir se existe alguma metodologia mais adequada para lidar com o tipo de incerteza do problema, se tal metodologia é computacionalmente tratável e quais as vantagens e desvantagens que as metodologias disponíveis na literatura podem trazer na análise do problema. Neste trabalho, são analisadas duas importantes metodologias para lidar com um problema de planejamento da produção sob incertezas: a programação estocástica de dois estágios e a otimização robusta. Ao passo que a programação estocástica é uma das técnicas mais tradicionalmente utilizadas em problemas de planejamento da produção sob incertezas, a mesma pode gerar modelos intratáveis se o número de cenários for muito grande. A otimização robusta surge como alternativa para superar a aparente limitação dos modelos de programação estocástica, mas ela pode ser muito conservadora, dependendo de como as incertezas são modeladas. As vantagens e desvantagens de cada metodologia são ilustradas com base num problema prático de planejamento da produção de empresas moveleiras e são comparadas em termos de função objetivo, nível de serviço e esforço computacional. Os resultados sugerem que ambas as técnicas são competitivas quando budgets de incerteza menos conservadores são utilizados no modelo de otimização robusta. Verificou-se também que o modelo equivalente robusto pode ser mais fácil de ser resolvido do que a versão estocástica, o que é especialmente importante quando a versão determinística já apresenta dificuldade de resolução.

Otimização sob incerteza; Planejamento da produção; Programação estocástica; Programação robusta

Optimizing production planning problems under uncertainty is a challenge, because it is necessary to define whether there is a methodology more appropriate to deal with the type of uncertainty of the problem, whether such methodology is computationally tractable, and which advantages and disadvantages the approaches available in the literature can bring to the analysis of the problem. In this paper, we analyze two important methodologies to deal with uncertainties in a production planning problem: two-stage stochastic programming and robust optimization. Whereas stochastic programming is one of the techniques most traditionally used in production planning problems under uncertainty, such approach can generate intractable models when a very large number of scenarios is considered. Robust optimization arises as an alternative technique to overcome the potential drawback of stochastic programming models, but it can be overly conservative depending on how the uncertainties are modeled. This paper also discusses the advantages and disadvantages of each methodology based on a practical problem of production planning in the furniture industry. The comparison between both approaches is given in terms of objective function, service level, and computational effort. The overall results suggest that both techniques are competitive when less conservative budgets of uncertainty are used in the robust optimization model. It was also verified that the robust counterpart model can be more easily solved compared with the stochastic version model, which is especially important when the deterministic model is already difficult to solve.

Optimization under uncertainty; Production planning; Stochastic programming; Robust optimization

1 Introdução

Na literatura de gestão de operações/pesquisa operacional, existem várias metodologias já consolidadas para lidar com problemas de otimização sob incertezas. Dentre as mais utilizadas, têm-se: análise de sensibilidade, programação estocástica e otimização robusta.

A análise de sensibilidade é um estudo de pós-otimalidade para determinar o impacto que as perturbações causam no problema nominal. Mulvey et al. (1995) consideram esse tipo de estudo reativo, uma vez que não é concebido para sugerir ações em uma situação de variação dos dados. A análise de sensibilidade pode não ser muito eficiente quando o interesse é analisar o impacto na solução ótima quando vários parâmetros variam simultaneamente. As outras metodologias mencionadas podem ser caracterizadas como pró-ativas no sentido de Mulvey et al. (1995). Basicamente, diferenciam-se pelos pressupostos teóricos e em como a incerteza é modelada (ou como ela é entendida em determinado contexto).

Na programação estocástica, por exemplo, assume-se que a distribuição de probabilidade dos parâmetros incertos (variáveis aleatórias) é conhecida ou pode ser razoavelmente bem estimada. O objetivo é determinar alguma política factível para todos – ou quase todos – os possíveis valores das realizações das variáveis aleatórias, de forma a otimizar algum funcional que depende das variáveis aleatórias. Uma importante restrição dessa metodologia é a suposição de que a distribuição de probabilidade dos dados é conhecida e não depende da decisão tomada (Kall & Wallace, 1994; Birge & Louveaux, 1997; Shapiro et al., 2009). Vale destacar que existem várias formulações possíveis dentro da programação estocástica, como modelos de dois estágios, modelos multiestágios, modelos com restrições probabilísticas (chance constraints), modelos estocásticos robustos e/ou com aversão ao risco, programação dinâmica estocástica, etc.

Finalmente, a otimização robusta preocupa-se em desenvolver modelos e métodos para que as soluções sejam factíveis para quaisquer realizações das variáveis aleatórias dentro de um conjunto convexo dado a priori e violações das restrições não são toleradas. Uma diferença bastante evidenciada entre essa metodologia e a programação estocástica é que a primeira não necessita do conhecimento da distribuição de probabilidade do parâmetro incerto, enquanto a segunda precisa (Ben-Tal & Nemirovski, 2000; Bertsimas & Sim, 2003).

Ao se optar por uma metodologia que lide com dados de entrada incertos, surgem algumas questões, como: “Existe alguma metodologia mais adequada para lidar com as incertezas de um determinado problema?”; “Tal metodologia é computacionalmente tratável para o problema que se quer resolver”?; ou “Quais as vantagens e desvantagens que as metodologias disponíveis na literatura podem trazer na análise do problema?”. Motivado por essas questões, este trabalho tem como objetivo analisar duas destas proeminentes metodologias para lidar com um problema de planejamento da produção sob incertezas: a programação estocástica de dois estágios e a otimização robusta. Além de ser a técnica mais utilizada em problemas de otimização sob incertezas, a programação estocástica é bastante versátil na incorporação de medidas de robustez/risco e na escolha dos estágios das variáveis de decisão, mas pode gerar modelos intratáveis se o número de cenários for muito grande. A otimização robusta, por sua vez, tem a vantagem de não utilizar cenários na descrição dos parâmetros incertos e gera modelos de mesma complexidade computacional que as versões determinísticas quando conjuntos de incerteza poliédricos são usados. Uma desvantagem da otimização robusta é a necessidade de garantir soluções ótimas (ou subótimas) para assegurar importantes propriedades teóricas. Para ilustrar a análise, as metodologias são aplicadas num problema prático de planejamento da produção de empresas moveleiras e são comparadas em termos de função objetivo, nível de serviço e esforço computacional.

O artigo está organizado da seguinte maneira. As Seções 2 e 3 apresentam, respectivamente, a fundamentação teórica dos métodos de programação estocástica e otimização robusta utilizadas no presente trabalho. A Seção 4 apresenta o problema de planejamento da produção de empresas moveleiras e desenvolve os modelos de programação estocástica e de otimização robusta. A Seção 5 ilustra e compara os resultados numéricos de ambas as abordagens. Finalmente, a Seção 6 descreve as principais conclusões.

2 Fundamentação teórica

2.1 Programação estocástica

Programação estocástica é uma das metodologias mais utilizadas para lidar com problemas de otimização sob incertezas. Na programação estocástica de dois estágios com recurso, tem-se um conjunto de variáveis de decisão de primeiro estágio que devem ser determinadas antes da realização das variáveis aleatórias. Por esse motivo, é comum atribuir a essas variáveis a denominação de decisões “aqui e agora” (here-and-now). Depois da realização das variáveis aleatórias, determinam-se as variáveis de decisão de segundo estágio, que são ações corretivas tomadas para adaptar ou finalizar as decisões feitas no primeiro estágio, de acordo com a observação dos valores das variáveis aleatórias. As variáveis de segundo estágio também são comumente chamadas de decisões “espere e veja” (wait-and-see) ou decisões de recurso (recourse decisions). O objetivo de um modelo de programação estocástica de dois estágios é identificar uma solução de primeiro estágio que seja bem equilibrada, diante de todas as possíveis realizações das variáveis aleatórias.

Em várias aplicações, é comum representar as variáveis aleatórias em algum espaço de probabilidade$(\Omega ,F,\Pi )$, em que$\Omega$ é o conjunto de possíveis estados da natureza (sendo que a realização genérica da variável aleatória $\xi$ é denotada por w) equipado com uma $\sigma -$álgebra de eventos F e com uma medida de probabilidade $\Pi$. Um modelo de programação estocástica de dois estágios com recurso ou, simplesmente, problema recurso RP, com incerteza apenas nos termos independentes (RHS) é apresentado na formulação matemática (1):

$$\begin{array}{l}\text{Minimizar}{c}^{T}x+\text{E}[q^{T}y(\omega )]\\ \text{Sujeito a:}Ax=b\\ Tx+Wy(\omega )\ge h(\omega )\\ x,y(\omega )\ge 0.\end{array}$$ (1)

No modelo (1), assume-se que os parâmetros c, q, A, b, T e W são parâmetros determinísticos e a função objetivo é caracterizada pelo valor esperado${\rm E}(\cdot )$. Para cada possível realização$\omega$, $h(\omega )$define o parâmetro estocástico referente ao RHS. Além disso, x é a variável de decisão de primeiro estágio e$y(\omega )$ define a variável de decisão de segundo estágio, como função da realização$\omega$.

É usual assumir que o parâmetro aleatório tem uma distribuição de probabilidade discreta com um número finito S de possíveis realizações${\omega }_s=[h_s]$, com suporte finito$\Omega =\{\mathrm{1,2,}\dots ,S\}$. Tais realizações são denominadas cenários e têm uma probabilidade de ocorrência que representa a chance de o evento materializar-se, dada por$\mathrm{Pr}(\omega ={\omega }_s)={\pi }_s$, de tal maneira que${\pi }_s>0$e${\displaystyle {\sum }_s{\pi }_s=1}$. Assim, o problema estocástico pode ser escrito na sua forma equivalente determinística, que resulta, em geral, num problema de otimização de grande porte, especialmente se o número de cenários for muito grande. Se para algum$s\in \Omega$, a restrição$Tx+W{y}_s\ge h_s$, com$x,y_s\ge 0$ não tiver solução, i.e., o correspondente problema de segundo estágio é infactível, então o equivalente determinístico do problema (1) é infactível. Nesse caso, assume-se que o valor da solução ótima é$+\infty$. O valor esperado de (1) é tomado em relação à distribuição de probabilidade da variável aleatória$\omega$, que é supostamente conhecida.

Mais recentemente, os problemas de programação estocástica como (1) têm sido estendidos para incorporar medidas de risco, cujo objetivo principal é controlar a variabilidade do custo total esperado de acordo com as preferências do decisor ou gerar soluções menos sensíveis às realizações dos cenários (ou mais “robustas”). Um dos primeiros trabalhos que se preocuparam com as questões de robustez em problemas de programação estocástica foi apresentado em Mulvey et al. (1995). Nesse trabalho pioneiro, os autores propõem uma nova metodologia baseada em programação estocástica (abordagem por cenários) e programação por metas, que foi designada pelos autores de otimização robusta. A motivação inicial dessa metodologia era desenvolver modelos de programação matemática cujas soluções permanecessem “próximas” da solução ótima e “quase” factíveis para quaisquer realizações das variáveis aleatórias. Soluções desses tipos seriam consideradas robustas em relação à otimalidade e à factibilidade, respectivamente.

Embora alguns autores ainda utilizem o termo otimização robusta ou programação estocástica robusta para designar tal metodologia, há uma tendência atual em utilizar o termo otimização aversa ao risco ou programação estocástica com aversão ao risco (risk-averse stochastic programming) para aludir à adoção de medidas de controle de robustez (ou medidas de risco) em modelos de programação estocástica (Alonso-Ayuso et al., 2014). No contexto dos problemas de programação estocástica com aversão ao risco, modelos do tipo média-risco tornaram-se populares, principalmente devido a Mulvey et al. (1995). A ideia principal desses modelos é combinar a expectância e a dispersão do custo total esperado, como ilustra a função objetivo (2):

Minimizar $\text{E}[g]+\lambda \cdot \text{D}[g]$(2)

em que g é o custo total do modelo de programação estocástica (1). A nova função objetivo (2) é caracterizada pelo valor esperado ${\rm E}(\cdot )$ e pela medida de dispersão $\text{D}(\cdot )$. O coeficiente $\lambda$ é um peso não negativo que controla o compromisso (tradeoff) entre custo esperado e risco (variabilidade). Aumentando-se $\lambda$, é possível gerar soluções com baixa variabilidade, mas com custos esperados elevados. O problema com $\lambda =0$representa o programa estocástico tradicional neutro ao risco (1). O leitor interessado em modelos média-risco pode consultar Ahmed (2006), que discute várias estratégias de minimização de risco.

Nos problemas de programação estocástica com aversão ao risco, é importante elaborar medidas de risco$\text{D}(\cdot )$ apropriadas. Há vários trabalhos que mencionam o modelo clássico de média-variância de Markowitz (Markowitz, 1959) a fim de balancear expectância e variabilidade, tais como Mulvey et al. (1995), Takriti & Ahmed (2004) e Zanjani et al. (2009). Medidas de risco baseadas em gestão financeira têm sido amplamente utilizadas na última década. Nesse contexto, vale destacar o valor no risco (VaR) e o condicional valor no risco (CVaR), introduzidos por Rockafellar & Uryasev (2000). Basicamente, o VaR representa a máxima perda esperada incorrida em um horizonte de planejamento a um nível de significância$\alpha$. O CVaR minimiza o valor no risco (VaR) juntamente com custo esperado excedendo VaR no nível de significância$\alpha$.

Embora não seja possível garantir a priori a melhor medida de risco possível de acordo com critérios pré-definidos, alguns autores têm proposto com sucesso modelos de programação estocástica com critérios de aversão ao risco para problemas específicos de planejamento da produção, como Alonso-Ayuso et al. (2014), Alem & Morabito (2013b), Toso & Alem (2014), Amorim et al. (2013), Barreiros et al. (2013), dentre outros.

2.2 Otimização robusta

A ideia da otimização robusta é desenvolver modelos “imunes” ou imunizados, tanto quanto for possível, às incertezas nos dados. Segundo esse viés, Soyster (1973) foi o trabalho pioneiro que propôs um modelo de otimização linear para gerar soluções factíveis para todo dado pertencente a um determinado conjunto convexo. O modelo de otimização resultante produz soluções ditas “muito conservadoras”, no sentido em que o valor da função objetivo deteriora-se em demasia para garantir a robustez – em termos de factibilidade – da solução. Sem perda de generalidade, considere o problema de otimização (3) com incertezas nos parâmetros do RHS do segundo conjunto de restrições:

$$\begin{array}{l}\text{Minimizar}{c}^{T}x+q^{T}y\\ \text{Sujeito a:}Ax=b\\ T^{a}x+Wy\ge h\\ x,y\ge 0.\end{array}$$ (3)

em que o segundo conjunto de restrição foi reescrito como $Tx-h{x}_{n+1}+Wy\ge 0$, com $x_{n+1}=1$ para incluir o parâmetro incerto h na matriz tecnológica T. Como resultado, temos a matriz tecnológica aumentada $T^a=[T|h]$.

Suponha que $J_i$ seja o conjunto dos coeficientes na linha i da matriz $T^a$ que estão sujeitos à incerteza. Cada entrada $t_{ij}^a$, $j\in J_i$ é modelada como uma variável aleatória simétrica e limitada ${\tilde{t}}_{ij}$, $j\in J_i$, que toma valores no intervalo $[t_{ij}^a-{\widehat{t}}_{ij},t_{ij}^a+{\widehat{t}}_{ij}]$. O parâmetro ${\widehat{t}}_{ij}$é denominado de desvio da variável aleatória. Associado ao dado incerto ${\tilde{t}}_{ij}$, define-se uma outra variável aleatória ${\eta }_{ij}=({\tilde{t}}_{ij}-t_{ij}^a)/{\widehat{t}}_{ij}$, que obedece a uma distribuição qualquer simétrica no intervalo $[-\mathrm{1,1}]$.

A formulação de Soyster é conservadora porque considera que a variável aleatória ${\tilde{t}}_{ij}$ é definida como ${\tilde{t}}_{ij}=t_{ij}^a+{\widehat{t}}_{ij}$ ou como ${\tilde{t}}_{ij}=t_{ij}^a-{\widehat{t}}_{ij}$, dependendo de qual situação representa o desvio mais desfavorável em relação às possíveis realizações da variável aleatória. Esse conservadorismo exacerbado motivou Bertsimas & Sim (2003, 2004) a proporem uma metodologia alternativa para controlar o grau de conservadorismo da solução robusta, pela introdução de um parâmetro que pode ser definido pelo usuário, por exemplo.

Considere a i-ésima restrição ${\displaystyle {\sum }_j{\tilde{t}}_{ij}{x}_j}+{\displaystyle {\sum }_j{w}_{ij}{y}_j}\ge h_i$. Para todo i, é introduzido o parâmetro ${\Gamma }_i$ (não necessariamente inteiro), denominado de budget de incerteza, tal que ${\Gamma }_i\in [\mathrm{0,}|J_i|]$. A função de ${\Gamma }_i$ é exatamente a de ajustar a robustez da solução em relação ao grau de conservadorismo. Intuitivamente falando, é improvável que todo ${\tilde{t}}_{ij}$ assuma o pior caso simultaneamente para todo $j\in J_i$. Então, o objetivo torna-se estar protegido de um número ${\Gamma }_i$ de coeficientes que o decisor define de acordo com as suas preferências e atitudes em relação ao risco, por exemplo. Considere a formulação não linear (4):

$$\begin{array}{l}\text{Minimizar}{c}^{T}x+q^{T}y\\ \text{Sujeito a:}Ax=b\\ {\displaystyle \sum _j{t}_{ij}^a{x}_j}+\underset{U}{\max }{\displaystyle \sum _{j\in S_i}{\widehat{t}}_{ij}{x}_j}+{\displaystyle \sum _j{w}_{ij}{y}_j}\ge h_i,\forall i\\ x,y\ge 0.\end{array}$$ (4)

em que $U=\left\{S_i\right|S_i\subseteq J_i,\left|S_i\right|\le {\Gamma }_i\}$.

Se ${\Gamma }_i=0$, então não há controle sobre as variáveis aleatórias e o problema resultante coincide com o problema determinístico original (nesse caso, note que a segunda parcela de (4) é zero); essa atitude é preferida por decisores neutros em relação ao risco. Se ${\Gamma }_i=\left|J_i\right|$, então o problema resultante coincide com o pior caso de Soyster, pois é permitido que todos os coeficientes assumam o pior caso simultaneamente. No meio termo, i.e., ${\Gamma }_i\in [\mathrm{0,}|J_i|]$, tem-se a flexibilidade de ajustar o tradeoffentre robustez e grau de conservadorismo da solução. Note que aumentar ${\Gamma }_i$ significa aumentar o conjunto de incerteza dado U.

Para reformular o modelo (4) como um modelo de otimização linear equivalente, o primeiro passo é definir a função de proteção da i-ésima restrição e mostrar que ela é equivalente à função objetivo de um problema de otimização linear. Para isso, considere que ${\Gamma }_i$ pode assumir apenas valores inteiros.

Proposição 1. Dado um vetor $x*$, a função de proteção da i-ésima restrição, representada pela equação (5),

$${\beta }_i(x*,{\Gamma }_i)=\underset{U}{\max }{\displaystyle \sum _{j\in S_i}{\widehat{t}}_{ij}}{x}_j^{\ast },$$ (5)

é igual à função objetivo do seguinte problema de otimização linear (6):

$$\begin{array}{l}{\beta }_i(x*,{\Gamma }_i)=\text{Maximizar}{\displaystyle \sum _{j\in J_i}{\widehat{t}}_{ij}}{x}_j^{\ast }{z}_{ij}\\ \text{Sujeito a:}{\displaystyle \sum _{j\in J_i}{z}_{ij}}\le {\Gamma }_i\\ \text{0}\le z_{ij}\le \mathrm{1,}\forall j\in J_i.\end{array}$$ (6)

Prova: Considerando que as únicas variáveis de decisão do problema (5) são $z_{ij},j\in J_i$, e que se trata de um problema de maximização, a solução será escolhida de forma a ter a maior quantidade possível de variáveis em 1. Assim, as ${\Gamma }_i$’s variáveis com maior coeficiente ${\widehat{t}}_{ij}$ serão 1 e as remanescentes serão zero.

Note que, sem o vetor-solução $x*$ em mãos, a função objetivo do problema (6) apresenta um produto entre duas variáveis de decisão. Entretanto, é possível “linearizá-la” a partir do correspondente problema dual. Para isso, considere os multiplicadores de Lagrange ${\theta }_i$ e ${\phi }_{ij}$ associados à primeira e segunda restrições do problema (6). Então, o problema dual (7) é definido da seguinte maneira:

$\begin{array}{l}\text{Minimizar}{\Gamma }_i{\theta }_i+{\displaystyle \sum _{j\in J_i}{\phi }_{ij}}\\ \text{Sujeito a:}{\theta }_i+{\phi }_{ij}\ge {\widehat{t}}_{ij}{x}_j^{\ast },\forall i,j\in J_i\\ {\theta }_i,{\phi }_{ij}\ge 0,\forall i,j\in J_i.\end{array}$(7)

Por dualidade forte, uma vez que o problema (6) é factível e limitado para todo ${\Gamma }_i\in [\mathrm{0,}|J_i|]$, então o problema dual (7) é também factível e limitado, e as suas funções objetivos coincidem. Usando a Proposição (2.1), tem-se que ${\beta }_i(x*,{\Gamma }_i)$ é igual ao valor da função objetivo do problema (7). Substituindo o vetor dado a priori $x*$ pela correspondente variável de decisão $x$ e injetando o dual (7) ao problema (4), obtém-se o equivalente robusto linear do problema nominal (3), que é representado no modelo (8):

$$\begin{array}{l}\text{Minimizar}{c}^{T}x+q^{T}y\\ \text{Sujeito a:}Ax=b\\ {\displaystyle \sum _j{t}_{ij}^a{x}_j}+{\theta }_i{\Gamma }_i+{\displaystyle \sum _{j\in J_i}{\phi }_{ij}}+{\displaystyle \sum _j{w}_{ij}{y}_j}\ge h_i,\forall i\\ x,y\ge 0\\ {\theta }_i+{\phi }_{ij}\ge {\widehat{t}}_{ij}{x}_j,\forall i,j\in J_i\\ {\theta }_i,{\phi }_{ij}\ge \mathrm{0,}\forall i,j\in J_i\end{array}$$ (8)

Quando ${\Gamma }_i$ for fracionário, a função de proteção da i-ésima restrição é definida de acordo com a equação (9):

$${\beta }_i(x*,{\Gamma }_i)=\underset{\{S_i\cup \{{\kappa }_i\left\}\right|S_i\subseteq J_i,\left|S_i\right|\le \lfloor {\Gamma }_i\rfloor ,{\kappa }_i\in J_i\backslash S_i\}}{\max }\left\{{\displaystyle \sum _{j\in S_i}{\widehat{t}}_{ij}}{x}_j^{\ast }+\left({\Gamma }_i-\lfloor {\Gamma }_i\rfloor \right){\widehat{t}}_{i{\kappa }_i}{x}_{{\kappa }_i}^{\ast }\right\},$$ (9)

e o desenvolvimento do equivalente robusto é análogo.

Por construção, o equivalente robusto é deterministicamente factível se, no máximo, ${\Gamma }_i$ coeficientes variarem. Porém, o que acontece se o decisor for “otimista” e escolher um valor baixo para ${\Gamma }_i$ (por exemplo, ${\Gamma }_i<<\left|J_i\right|$), mas na prática ele observar que mais coeficientes variaram (assumiram) valores desfavoráveis? Nesse caso, a metodologia de otimização robusta fornece alguns limitantes para determinar a probabilidade de violação do primeiro conjunto de restrição em (8). Seja $x*$ a solução do modelo de otimização robusta. Então, a probabilidade de violação da i-ésima restrição pode ser aproximada pela inequação (10):

$$\mathrm{Pr}\left({\displaystyle \sum _j{\tilde{t}}_{ij}{x}_j^{\ast }}+{\displaystyle \sum _y{w}_{ij}{y}_j}<h_i\right)\le 1-\Phi \left(\frac{{\Gamma }_i-1}{\sqrt{\left|J_i\right|}}\right)$$,(10)

em que $\Phi (\eta )=\frac{1}{\sqrt{2\pi }}{\displaystyle \underset{-\infty }{\overset{\eta }{\int }}\exp -\left(\frac{v^2}{2}\right)}dv$ é a função de distribuição acumulada de uma normal padrão.

A partir da expressão (8) é possível determinar valores para os budgets de incerteza. Note que se $\mathrm{Pr}({\displaystyle {\sum }_j{\tilde{t}}_{ij}{x}_j^{\ast }+}{\displaystyle {\sum }_j{w}_{ij}{y}_j}<h_i)={\epsilon }_i$, então a inequação (10) pode ser escrita como mostra a equação (11).

$${\epsilon }_i=1-\Phi \left(\frac{{\Gamma }_i-1}{\sqrt{\left|J_i\right|}}\right)\Rightarrow \frac{{\Gamma }_i-1}{\sqrt{\left|J_i\right|}}={\Phi }^{-1}(1-{\epsilon }_i)\Rightarrow {\Gamma }_i=1+{\Phi }^{-1}(1-{\epsilon }_i)\sqrt{\left|J_i\right|}$$ (11)

O leitor interessado em aplicações de otimização robusta em planejamento e programação da produção pode consultar Ben-Tal et al. (2005); Bertsimas & Thiele (2006); Li & Ierapetritou (2008); Palma & Nelson (2009); Bohle et al. (2009); Alem & Morabito (2012, 2013a); Munhoz & Morabito (2013); Paiva & Morabito (2014), dentre outros.

3 Aplicação: planejamento da produção em fábricas de móveis

3.1 Modelo determinístico

Nessa seção, é apresentado um modelo matemático inteiro-misto para o problema de planejamento da produção típico de empresas moveleiras (Gramani & França, 2006; Rangel & Figueiredo, 2008; Alem & Morabito, 2012). A ideia básica do modelo é combinar as decisões de dimensionamento dos lotes de produção – quais móveis produzir, quanto e quando produzir – às decisões operacionais de corte de estoque, e.g., quais padrões de corte serão usados e a frequência de utilização para produzir as peças necessárias na montagem dos móveis. Sejam os conjuntos I de produtos (guarda-roupas, cômodas, criados-mudos, camas, etc.) e P das peças que compõem os produtos que devem ser produzidas de acordo com um conjunto J de padrões de corte, de modo a atender à demanda ao longo de um conjunto T de períodos do horizonte de planejamento. A Tabela 1 lista a notação usada na formulação matemática do problema.

Tabela 1
Notação matemática dos modelos de otimização. Fonte: Elaborada pelos autores.

Dados de Entrada
$c_{it}$	Custo de produção do produto i no período t
$h_{it}^{+}$	Custo de estoque do produto i no período t
$h_{it}^{-}$	Custo de atraso do produto i no período t
$w_{jt}$	Custo de desperdício de material ao se cortar o padrão de corte j no período t
$c{s}_{it}$	Custo de preparar a máquina para produzir o produto i no período t
$o_t$	Custo de hora extra no período t
$a_{pj}$	Número de vezes que a peça p aparece no padrão de corte j
$r_{pi}$	Número de peças p que compõem o produto i
$d_{it}$	Demanda do produto i no período t
$C_t$	Capacidade regular no período t
$C_t^E$	Capacidade extra no período t
$I_{it}^{\max }$	Estoque máximo permitido do produto i no período t
$v_i$	Tempo de produção do produto i
$M$	Número suficientemente grande
Variáveis de decisão
$O_t$	Hora extra utilizada no período t
$X_{it}$	Quantidade de produto i produzido no período t
$I_{it}^{+}$	Quantidade de produto i em estoque ao final do período t
$I_{it}^{-}$	Quantidade de produto i em atraso ao final do período t
$Y_{jt}$	Frequência do padrão de corte j no período t
$Z_{it}$	Variável binária que vale 1 se produto i for produzido no período t

O modelo matemático inteiro-misto a seguir objetiva determinar um plano de produção de móveis a um custo mínimo, a partir do processamento de padrões de corte pré-selecionados, de modo a atender à demanda e satisfazer as restrições do processo produtivo. A seguir, a primeira formulação matemática do problema combinado de dimensionamento de lotes e corte de estoque CLC é detalhada.

A função objetivo (12) consiste em minimizar o custo total de produção, estoque, atraso, preparação, perda de material e hora extra. O primeiro termo em (12) é o custo total incorrido na produção, estoque, atraso e preparação. O segundo termo representa o custo de desperdício de material. O último termo é o custo devido à utilização de horas extras.

$$\text{Minimizar}{\displaystyle \sum _{i\in I}{\displaystyle \sum _{t\in T}\left(c_{it}{X}_{it}+h_{it}^{+}{I}_{it}^{+}+h_{it}^{-}{I}_{it}^{-}+c{s}_{it}{Z}_{it}\right)}}+{\displaystyle \sum _{j\in I}{\displaystyle \sum _{t\in T}{w}_{jt}{Y}_{it}}}+{\displaystyle \sum _{t\in T}{o}_t{O}_t}$$ (12)

As restrições de balanceamento de estoque (13)fazem o balanço de toda a produção de móveis. Sem perda de generalidade, assume-se que os níveis de estoque e atraso no início do horizonte de planejamento são nulos.

$$X_{it}+I_{it-1}^{+}+I_{it}^{-}+=d_{it}+I_{it}^{+}+I_{it-1}^{-},\forall i\in I,t\in T.$$ (13)

As restrições de balanceamento de peças (14)forçam o balanço de peças necessárias para montar os produtos. Note que essas restrições de acoplamento são as únicas que integram ambas as decisões relacionadas ao dimensionamento de lotes $X_{it}$ e ao corte de estoque $Y_{jt}$.

$${\displaystyle \sum _{j\in J}{a}_{pj}}{Y}_{jt}\ge {\displaystyle \sum _{i\in I}{r}_{pi}}{X}_{it},\forall p\in P,t\in T.$$ (14)

As restrições de estoque (15) limitam o máximo estoque de móveis permitido e o mínimo estoque de móveis requerido. Nesse caso, o mínimo estoque é zero, mas quantidades positivas podem ser consideradas, dependendo da política da empresa (por exemplo, estoque de segurança).

$$0\le I_{it}^{+}\le I_{it}^{\max },\forall i\in I,t\in T.$$ (15)

As restrições de capacidade (16) indicam que o tempo total de produção dos móveis deve ser inferior à soma das capacidades regular e extra. Pode-se também adotar outra unidade de capacidade de produção, como a área total cortada em metros quadrados, por exemplo.

$${\displaystyle \sum _{i\in I}{v}_i}{X}_{it}\le C_t+O_t,\forall t\in T.$$ (16)

As restrições de hora extra (17) limitam a utilização da quantidade de horas extras em cada período. Pode-se considerar a capacidade extra como uma fração da capacidade regular.

$$0\le O_t\le C_t^E,\forall t\in T.$$ (17)

As restrições de preparação (18) indicam que pode haver produção do produto i no período t somente se a linha de produção estiver preparada.

$$X_{it}\le M{Z}_{it},\forall i\in I,t\in T.$$ (18)

O conjunto de restrições (19)-(22) refere-se ao domínio das variáveis de decisão.

$$X_{it}\ge \mathrm{0,}\forall i\in I,t\in T.$$ (19)

$$I_{it}^{-}\ge \mathrm{0,}\forall i\in I,t\in T.$$ (20)

$$Y_{jt}\ge \mathrm{0,}\forall j\in J,t\in T.$$ (21)

$$Z_{it}\in \left\{\mathrm{0,1}\right\},\forall i\in I,t\in T.$$ (22)

Neste trabalho, as variáveis relacionadas aos padrões de corte, $Y_{jt}$, são consideradas números não negativos reais, mas é possível impor restrições de integralidade sobre elas para se obter apenas quantidades inteiras. Vários trabalhos da literatura excluem a restrição de integralidade sobre tais variáveis na tentativa de obter um modelo de otimização mais tratável computacionalmente. Uma justificativa é que o simples arredondamento destas variáveis na solução do problema é razoável em vários contextos práticos.

3.2 Modelo de programação estocástica

Nessa seção, é desenvolvido um modelo de planejamento da produção sob incerteza em fábricas de móveis sob o viés da programação estocástica de dois estágios com recurso em que a demanda éstocástica. Note que, nesse caso, a demanda estocástica corresponde ao parâmetro incerto RHS do modelo 1 apresentado na Seção 2.1. Para determinar o modelo determinístico equivalente, assume-se que a demanda estocástica é aproximada por um conjunto discreto e finito de cenários, como já discutido. Na formulação de programação estocástica, o estoque $I_{its}^{+}$ e o atraso $I_{its}^{-}$ foram considerados variáveis de decisão de segundo estágio, sendo ponderadas na função objetivo pela probabilidade ${\pi }_s$ associada a cada cenário. As outras variáveis de decisão foram designadas como primeiro estágio. Além disso, introduziu-se um parâmetro de risco ${\Delta }_s$ para controlar o desvio do custo total esperado em relação ao custo de segundo estágio de cada cenário $s\in \Omega$. A versão estocástica com aversão ao risco do tipo média-risco do modelo CLC é chamada, a partir desse ponto, de modelo MR, e pode ser escrita como mostra a função objetivo (23):

$$\text{Minimizar}{\displaystyle \sum _{i\in I}{\displaystyle \sum _{t\in T}\left(c_{it}{X}_{it}+c{s}_{it}{Z}_{it}\right)}}+{\displaystyle \sum _{j\in I}{\displaystyle \sum _{t\in T}{w}_{jt}{Y}_{it}}}+{\displaystyle \sum _{i\in I}{\displaystyle \sum _{t\in T}{\displaystyle \sum _{s\in \Omega }{\pi }_s\left(h_{it}^{+}{I}_{its}^{+}+h_{it}^{-}{I}_{its}^{-}\right)+}}}{\displaystyle \sum _{t\in T}{o}_t{O}_t}+\lambda {\displaystyle \sum _{s\in \Omega }{\pi }_s}{\Delta }_s$$ (23)

Sujeito a: Restrições (14), (16)-(19), (21), (22)

$$\begin{array}{l}{X}_{it}+I_{i(t-1)s}^{+}+I_{its}^{-}+=d_{its}+I_{its}^{+}+I_{i(t-1)s}^{-},\forall i\in I,t\in T,s\in \Omega .\\ 0\le I_{its}^{+}\le I_{it}^{\max },\forall i\in I,t\in T,s\in \Omega \\ {\Delta }_s\ge {\displaystyle \sum _{i\in I}{\displaystyle \sum _{t\in T}\left(h_{it}^{+}{I}_{its}^{+}+h_{it}^{-}{I}_{its}^{-}\right)-}}{\displaystyle \sum _{i\in I}{\displaystyle \sum _{t\in T}{\displaystyle \sum _{s\text{'}\in \Omega }{\pi }_{s\text{'}}\left(h_{it}^{+}{I}_{its\text{'}}^{+}+h_{it}^{-}{I}_{its\text{'}}^{-}\right)}}},s\in \Omega \\ I_{its}^{-}\ge \mathrm{0,}\forall i\in I,t\in T,s\in \Omega \\ {\Delta }_s\ge \mathrm{0,}s\in \Omega \end{array}$$ (24)

As restrições (24) são similares às suas versões determinísticas, mas devem ser satisfeitas para todo cenário $s\in \Omega$. O modelo também tem um caráter multicritério com dois objetivos explícitos na função objetivo: a primeira parcela representa o custo total médio ou esperado (soma dos custos de primeiro e segundo estágios) que ignora o custo de variabilidade, enquanto a segunda parcela é exatamente o custo da variabilidade ponderado pelo parâmetro de risco $\lambda$. Quanto maior for o desvio ${\Delta }_s$, maior é o risco atribuído à solução para um dado s, uma vez que, nesse caso, a solução do cenário s é superior ao custo total esperado. Então, ao se minimizar o desvio, os custos dos diferentes cenários tendem a se aproximar e, com isso, a solução ótima de um cenário spermanece “perto” da solução ótima para qualquer outro cenário s’, e a solução pode ser dita robusta (Mulvey et al., 1995). Para $\lambda =0$, a parcela referente ao risco é ignorada. A princípio, não há um valor máximo para $\lambda$. Diz-se, simplesmente, que quanto maior é esse valor, mais o decisor é averso ao risco, no sentido de que ele está disposto a sacrificar o custo total esperado em troca obter soluções mais robustas. Em geral, são resolvidos vários modelos MR variando-se $\lambda$ num certo intervalo na tentativa de determinar as curvas de tradeoff entre custo total e robustez. Assim, o decisor pode fazer a escolha de $\lambda$ de acordo com a sua atitude em relação ao risco.

3.3 Modelo de otimização robusta

Para aplicar a metodologia de otimização robusta ao CLC com demanda incerta, reescreve-se a restrição (13), pois ela apresenta uma igualdade, que pode não ser satisfeita para todas as realizações da demanda. Assim, a restrição (13) é redefinida como uma restrição de custo de estoque e atraso, da seguinte maneira. Defina o estoque do produto i no período t como $I_{it}=I_{it}^{+}-I_{it}^{-}$. Note que o estoque definido dessa maneira é irrestrito em sinal, podendo representar estoque em mãos ou atraso. Uma vez que $I_{it}={\displaystyle \sum _{\tau =1}^t{X}_{i\tau }}-d_{i\tau }$ e definindo as variáveis não negativas de custo ou atraso, $H_{it}$, obtém-se as inequações (25) e (26):

$$H_{it}\ge h_{it}^{+}{I}_{it}=h_{it}^{+}\left({\displaystyle \sum _{\tau =1}^t(X_{i\tau }}-d_{i\tau })\right),\forall i\in I,t\in T$$ (25)

$$H_{it}\ge h_{it}^{-}(-I_{it})=h_{it}^{-}\left(-{\displaystyle \sum _{\tau =1}^t(X_{i\tau }}-d_{i\tau })\right),\forall i\in I,t\in T$$ (26)

Agora, pode-se aplicar as técnicas de otimização robusta para lidar com a demanda incerta nas restrições reformuladas (25) e (26). Assim, considere que as demandas acumuladas sejam variáveis aleatórias simétricas e limitadas ${\tilde{d}}_{i\tau }$ que assumem valores no intervalo $\left[d_{i\tau }-{\widehat{d}}_{i\tau },d_{i\tau }+{\widehat{d}}_{i\tau }\right]$. O desvio relativo entre as demandas incertas e nominais é determinado por $z_{i\tau }=({\tilde{d}}_{i\tau }-d_{i\tau })/{\widehat{d}}_{i\tau }$ que pertence ao intervalo $[-\mathrm{1,1}]$. Baseado no trabalho de Bertsimas & Thiele (2006), os budgets de incerteza${\Gamma }_{it}\in [\mathrm{0,}t]$ são assumidos não decrescentes com o período t. Além disso, supõe-se que o aumento em ${\Gamma }_{it}$ é igual ou menor que o aumento nos períodos.

De uma perspectiva de pior caso, deve-se maximizar o lado direito das restrições (25) e (26) sobre o conjunto de todas as possíveis realizações das demandas incertas. Para i e t dados, isso corresponde a resolver o problema auxiliar (27):

$$\begin{array}{l}\text{Maximizar}{\displaystyle \sum _{\tau =1}^t{\widehat{d}}_{i\tau }}{z}_{i\tau }\\ \text{Sujeito a:}{\displaystyle \sum _{\tau =1}^t{z}_{i\tau }}\le {\Gamma }_{it}\\ \text{0}\le z_{i\tau }\le \mathrm{1,}\forall \tau \le t\end{array}$$ (27)

De acordo com a filosofia da otimização robusta, o problema auxiliar (27) resulta em minimizar${\displaystyle \sum _{\tau =1}^t{\widehat{d}}_{i\tau }}{z}_{i\tau }$ na restrição (25) e maximizar ${\displaystyle \sum _{\tau =1}^t{\widehat{d}}_{i\tau }}{z}_{i\tau }$na restrição (26). Similarmente à Seção 2.2, pode-se determinar o equivalente robusto do problema (12)-(22) com demanda incerta (denominado de OR), que é representado pelo modelo (28):

$$\text{Minimizar}{\displaystyle \sum _{i\in I}{\displaystyle \sum _{t\in T}\left(c_{it}{X}_{it}+h_{it}^{+}{I}_{it}^{+}+h_{it}^{-}{I}_{it}^{-}+c{s}_{it}{Z}_{it}\right)}}+{\displaystyle \sum _{j\in I}{\displaystyle \sum _{t\in T}{w}_{jt}{Y}_{it}}}+{\displaystyle \sum _{t\in T}{o}_t{O}_t}$$

Sujeito a: Restrições (14), (16)-(19), (21), (22)

$$\begin{array}{l}{H}_{it}\ge h_{it}^{+}\left({\displaystyle \sum _{\tau =1}^t(X_{i\tau }}-d_{i\tau })+{\theta }_{it}{\Gamma }_{it}+{\displaystyle \sum _{\tau =1}^t{\phi }_{i\tau t}}\right),\forall i\in I,t\in T\\ H_{it}\ge h_{it}^{-}\left({\displaystyle \sum _{\tau =1}^t-(X_{i\tau }}-d_{i\tau })+{\theta }_{it}{\Gamma }_{it}+{\displaystyle \sum _{\tau =1}^t{\phi }_{i\tau t}}\right),\forall i\in I,t\in T\\ {\displaystyle \sum _{\tau =1}^t(X_{i\tau }}-d_{i\tau })+{\theta }_{it}{\Gamma }_{it}+{\displaystyle \sum _{\tau =1}^t{\phi }_{i\tau t}}\le I_{it}^{\max },\forall i\in I,t\in T\\ {\theta }_{it}+{\phi }_{i\tau t}\ge {\widehat{d}}_{it}^{},\forall i\in I,t\in T,\tau \le t\\ {\theta }_{it},{\phi }_{i\tau t}\ge \mathrm{0,}\forall i\in I,t\in T,\tau \le t\end{array}$$ (28)

em que ${\theta }_{it}$ e ${\phi }_{i\tau t}$ são as variáveis duais associadas ao dual da formulação (27). A formulação (28) assegura que $H_{it}\ge 0$, pois se trata do custo de estoque e de falta, que são sempre não negativos. No trabalho de Wei et al. (2011), os autores derivam limitantes probabilísticos de violação de restrições de estoque e de falta. Considerando que $X_{it}^{\ast }$ seja a solução do equivalente robusto com demanda incerta, tem-se que a probabilidade de violação das restrições podem ser escritas como as inequações (29) e (30):

$$\mathrm{Pr}\left[H_{it}<h_{it}^{+}\left({\displaystyle \sum _{\tau =1}^t(X_{it}^{\ast }-{\tilde{d}}_{i\tau })}\right)\right]\le 1-\Phi \left(\frac{{\Gamma }_{it}-1}{\sqrt{t}}\right)$$,(29)

e

$$\mathrm{Pr}\left[H_{it}<h_{it}^{-}\left({\displaystyle \sum _{\tau =1}^t-(X_{it}^{\ast }-{\tilde{d}}_{i\tau })}\right)\right]\le 1-\Phi \left(\frac{{\Gamma }_{it}-1}{\sqrt{t}}\right)$$,(30)

em que $\Phi (\cdot )$ é a função de distribuição acumulada de uma normal padrão, para todo i, t.

4 Exemplo numérico ilustrativo

Nessa seção, são apresentados os resultados referentes aos dois modelos de otimização sob incertezas propostos. Os modelos foram codificados no Sistema de Modelagem Algébrica GAMS (McCarl et al., 2010) e resolvidos softwareILOG-CPLEX 11.0 (ILOG, 2008), usando-se os valores default para todos os parâmetros do CPLEX. Os experimentos foram executados num notebook Core 2 Duo 4, 4.0 GB de memória RAM, 2.0 GHz, sob a plataforma Windows VISTA. Os dados de entrada utilizados referem-se a informações reais de uma planta moveleira de pequeno porte estudada em Alem & Morabito (2012, 2013a, b).

Foram consideradas 3 famílias de produtos: guarda-roupas de 5 portas (a5p), cômodas (cmd) e criados-mudos (crd), os quais são produzidos a partir do corte de placas de MDF que são adquiridas dos varejistas em duas dimensões 2,75 x 1.83 m e 2,75 x 1,85 m e em 6 espessuras distintas: 3, 9, 12, 15, 20 e 25 mm. Os produtos são formados por um total de 49 peças retangulares. Tem-se um total de 81 padrões de corte, que são frequentemente utilizados pela empresa. A demanda refere-se a um período de 8 semanas de produção. Os detalhes relativos aos dados de entrada podem ser consultados em Alem (2011).

Para o problema de programação estocástica MR, foram gerados 100 cenários, segundo uma distribuição uniforme discreta entre $[\mathrm{0,7}{d}_{it};\mathrm{1,3}{d}_{it}]$ (para considerar um nível de variabilidade 30% abaixo e acima da demanda nominal) e todos foram admitidos equiprováveis. Para o modelo de otimização robusta OR, adotou-se ${\Gamma }_{it}=g\cdot t$, com $g=\mathrm{0,01}$ até 1. Dessa forma, vários budgets de incerteza foram testados, desde os menos conservadores até o mais conservador (${\Gamma }_{it}=t$). Na comparação, considerou-se que o budget mais conservador implica a solução mais robusta. Assim, quando ${\Gamma }_{it}=t$, o problema é totalmente robusto ou 100%; quando ${\Gamma }_{it}=\mathrm{0,5}t$, o problema é 50% robusto, e assim por diante. Para o modelo MR, a robustez do problema é representada pela porcentagem de decréscimo do desvio esperado (${\displaystyle {\sum }_s{\pi }_s{\Delta }_s}$) em relação ao desvio para $\lambda =0$ (problema estocástico tradicional). Quando o desvio atinge 30%, por exemplo, diz-se que o problema é 30% robusto.

Os modelos foram analisados de acordo com: (1) o aumento no valor ótimo da função objetivo; (2) o nível de serviço tipo II ou taxa de atendimento da demanda; e (3) o tempo de execução computacional. Tais critérios de comparação são baseados em trabalhos da literatura. De fato, analisar a deterioração no valor da função objetivo é uma prática comum dos trabalhos que desenvolvem modelos de otimização robusta e/ou modelos de programação estocástica com medidas de aversão ao risco (Bertsimas & Sim, 2004; Alonso-Ayuso et al., 2014, entre muitos outros). A motivação é analisar se a redução do risco (chance de violação das restrições e/ou variabilidade dos custos) “compensa” o custo extra incorrido na função objetivo. Da mesma forma, a análise do tempo computacional faz-se importante porque os modelos de programação estocástica têm a fama de se tornarem intratáveis quando o número de cenários aumenta demasiadamente. Finalmente, o nível de serviço foi analisado por ser, de certa forma, uma medida que reflete a qualidade da solução dos planos de produção obtidos em cada um dos modelos. De fato, soluções com baixos níveis de serviço são ruins pois implicam que parte significativa da demanda não pode ser atendida.

Para o problema de programação estocástica, o nível de serviço em cada cenário s foi determinado como: ${{\rm B}}_s=1-I_s^{-}/D_s$, em que $I_s^{-}={\displaystyle {\sum }_{i,t}{I}_{its}^{-}}$. A demanda total $D_s$ foi computada como $D_s={\displaystyle {\sum }_{i,t}{d}_{its}}$. O nível de serviço global do problema foi calculado como ${\displaystyle {\sum }_s{\pi }_s{{\rm B}}_s}$. Como o interesse foi determinar apenas a fração da demanda não atendida ao final do horizonte de planejamento, analisou-se apenas $I_{iTs}^{-}$, em que T é o último período do horizonte de planejamento. Para o problema de otimização robusta, $I_{iT}^{-}$ foi determinado com base no custo de atraso $H_{iT}$ e o cálculo do nível de serviço foi calculado de maneira análoga.

4.1 Resultados e discussão

Aumento no valor ótimo. As Figuras 1 e 2 ilustram o aumento no valor ótimo da função objetivo e a robustez dos modelos MR e OR, de acordo com o aumento do fator de risco ($\lambda$) e dos budgets de incerteza (${\Gamma }_{it}$), respectivamente. Da Figura 1, tem-se que o aumento no valor ótimo é mais acentuado até $\lambda =10$, decrescendo entre (10,33) e atingindo o menor crescimento a partir de $\lambda =33$. Em $\lambda =10$, observa-se que o modelo é 60% robusto, de acordo com a definição proposta. A partir desse ponto, a robustez aumenta lentamente numa transição de fases. Para alcançar um modelo totalmente robusto, incorre-se num aumento do valor ótimo de 120%. Por outro lado, na Figura 2, o aumento no valor ótimo é bastante lento até, aproximadamente, ${\Gamma }_{it}=\mathrm{0,15}t$. Nesse ponto, atinge-se 15% de robustez a uma deterioração de 20% do valor ótimo. Os próximos budgets, entretanto, são bastante conservadores; por exemplo, para atingir uma robustez de 50%, o valor ótimo aumenta 375%. O aumento no valor ótimo segue de forma significativa até ${\Gamma }_{it}=\mathrm{0,7}t$, a partir do qual é menos acentuado até $g=\mathrm{0,95}$, voltando a se acentuar no pior caso. Nesse cenário pessimista, 100% de robustez é alcançada por um aumento de quase 800%.

Figura 1
Aumento no valor ótimo da função objetivo e robustez do modelo MR, considerando o fator de risco entre 0 e 90. Fonte: Elaborada pelos autores.

Figura 2
Aumento no valor ótimo e robustez do modelo OR, considerando os budgets de incerteza ${\Gamma }_{it}=g\cdot t$, com g = 0,01 até 1. Fonte: Elaborada pelos autores.

Assim, pode-se inferir que, em relação ao aumento no valor ótimo, o modelo OR é competitivo com o modelo MR apenas para budgets de incerteza pouco conservadores. Por exemplo, no modelo MR, 11% de robustez é atingida a um custo 32% maior; no modelo OR a mesma robustez resulta em um aumento do valor ótimo de apenas 11%. Situação análoga pode ser verificada para alcançar uma robustez de 15%. Para budgets de incerteza intermediários e conservadores, o modelo OR é superado pelo modelo MR, uma vez que MR consegue melhores valores de robustez com custos bem menores, i.e., possui melhores tradeoffs entre custo e robustez. Esses resultados confirmam os testes anteriores sobre o desempenho pobre do modelo de otimização robusta OR com budgets pessimistas, ocasionado pela “superproteção” da restrição de atendimento da demanda Alem & Morabito (2012, 2013a).

Nível de serviço. As Figuras 3 e 4 exibem o comportamento do nível de serviço à medida que a robustez aumenta, considerando: (a) apenas a demanda perdida (i.e. demanda atrasada no último período do horizonte de planejamento); e (b) a soma da demanda atrasada até o último período (composta pela demanda perdida em T e a demanda atrasada em $\forall t=\mathrm{1..}T-1$. Em relação ao nível de serviço, note que o modelo OR tem vantagem sobre o modelo estocástico até, aproximadamente, g = 0,15; nesses casos, o nível de serviço do modelo de otimização robusta é 100%. Todavia, o modelo MR mantém ótimos níveis de serviço (maiores do que 93%) para todos os valores de robustez, ao passo que o nível de serviço do modelo OR deteriora-se razoavelmente quando a solução se torna mais robusta. Analisando o nível de serviço referente ao total de demanda atrasada e perdida, a situação é ainda mais vantajosa para o modelo MR, que ainda mantém boas taxas de atendimento da demanda (cerca de 80%) – diferentemente do modelo OR – cujo nível de serviço cai significativamente quando as demandas atrasadas são contabilizadas.

Figura 3
Nível de serviço do modelo MR de acordo com o aumento da robustez. Fonte: Elaborada pelos autores.

Figura 4
Nível de serviço do modelo OR de acordo com o aumento da robustez. Fonte: Elaborada pelos autores.

Tempo de execução. As Figuras 5 e 6 ilustram os tempos de execução do algoritmo branch-and-cut do CPLEX 11.0 para resolver até a prova de otimalidade os exemplares do modelo MR e OR, respectivamente. Como já era de se esperar, o modelo MR é mais difícil de ser resolvido do que o modelo OR. Basicamente, tem-se duas razões que podem concorrer para aumentar a dificuldade do modelo MR1: o número de cenários e a parcela referente ao risco para $\lambda >0$ (para $\lambda =0$, tem-se o modelo estocástico tradicional. De fato, a Tabela 2 mostra que o exemplar resolvido do modelo MR tem mais restrições e variáveis do que o modelo OR, o que compromete o seu tempo de execução. Além disso, o tempo de execução do modelo MR parece aumentar, à medida que mais robustez é imposta ao modelo, o que não parece acontecer ao modelo OR. Na Figura 5, é possível observar que o tempo médio de resolução do modelo MR foi 47 segundos, contra 0,36 segundos da versão OR. Considerando o exemplar resolvido, o modelo OR é superior ao modelo MR, em termos de tempo de execução. Embora o exemplar analisado seja pequeno, a potencial intratabilidade dos modelos de programação estocástica com o crescimento do número de cenários é bem conhecida na literatura.

Figura 5
Tempo de execução computacional do modelo MR de acordo com o aumento da robustez. Fonte: Elaborada pelos autores.

Figura 6
Tempo de execução computacional do modelo OR de acordo com o aumento da robustez. Fonte: Elaborada pelos autores.

Tabela 2
Número de variáveis de decisão e restrições de um exemplar dos modelos MR e OR. Fonte: Elaborada pelos autores.

Modelo	Restrições	Variáveis contínuas	Variáveis binárias
MR	4516	6372	24
OR	688	896	24

Assim, a partir dessa análise, é possível afirmar que o modelo MR e o modelo OR com budgets de incerteza pouco conservadores podem ser usados por tomadores de decisão aversos ao risco, pois ambos garantem a robustez do modelo a um custo relativamente baixo. Entretanto, se o nível de robustez desejado for muito elevado, budgetsde incerteza mais conservadores precisam ser adotados, e o modelo OR pode ser desvantajoso em relação ao modelo MR. Nesse caso, pode-se utilizar o modelo MR com um fator de risco elevado.

De forma geral, ambas as metodologias são valiosas para lidar com problemas de planejamento da produção sob incertezas. Cada uma apresenta diferenças marcantes em relação à maneira como a incerteza é representada, assim como suposições, simplificações e limitações próprias. O método mais adequado depende da aplicação específica e do tipo de dados que o decisor tem em mãos. Considerando a aplicação proposta neste trabalho, é possível citar algumas vantagens e limitações de cada metodologia, que podem ser úteis também na escolha do método para tratar problemas similares e correlatos.

As principais vantagens e potenciais desvantagens/limitações dos modelos de programação estocástica identificadas neste trabalho são:

1. i

   Permitem modelar as variáveis aleatórias de forma bastante natural com a utilização de cenários, os quais podem ser gerados com diferentes estruturas e distribuições de probabilidade. Para o tomador de decisão, é importante ter a possibilidade de incorporar ou desconsiderar cenários com padrões totalmente distintos, uma vez que essa estratégia pode flexibilizar o plano de produção.

2. ii

   A possibilidade de utilizar as variáveis de decisão de segundo estágio é interessante porque tais variáveis podem ser usadas para se proteger dos parâmetros estocásticos, ajustando e corrigindo decisões de primeiro estágio tomadas antes da realização das variáveis aleatórias. Além disso, é possível testar várias configurações para os conjuntos de variáveis de primeiro e segundo estágios.

3. iii

   O modelo de média-risco propicia uma análise multiobjetivo relativamente simples entre o custo da função objetivo e a robustez do modelo, uma vez que basta designar valores para o fator de risco e analisar as soluções geradas. Dessa forma, os tomadores de decisão podem incorporar ao modelo a sua atitude em relação ao risco.

4. iv

   Como não existe uma fórmula geral para dividir as variáveis de primeiro e segundo estágios, cabe ao decisor escolher quais vão pertencer a um ou a outro grupo, o que pode provocar confusão e gerar modelos incoerentes.

5. v

   (a) Assumir uma distribuição de probabilidade para as variáveis aleatórias; (b) atribuir as probabilidades dos diferentes cenários; e (c) decidir sobre a estrutura da árvore de cenários e/ou a quantidade de cenários que devem ser considerados, são tarefas em geral difíceis de serem realizadas na prática, e requerem conhecimento profundo do problema a ser tratado e de técnicas especializadas de resolução.

6. vi

   Dependendo da escolha da estrutura da árvore de cenários e da quantidade de cenários necessários para representar as variáveis aleatórias, o modelo estocástico pode se tornar computacionalmente intratável.

As principais vantagens e potenciais desvantagens/limitações dos modelos de otimização robusta identificadas neste trabalho são:

1. i

   Não é necessário inferir distribuições de probabilidade para modelar os parâmetros incertos, uma vez que os equivalentes robustos só requerem que as variáveis aleatórias sejam simétricas num intervalo pré-definido. Para o tomador de decisão, pode ser bastante natural inferir um limitante mínimo e máximo para a variação dos parâmetros incertos, baseados em dados históricos ou mesmo na sua experiência, em vez de determinar exatamente o valor da variável em cada cenário e a sua probabilidade de ocorrência.

2. ii

   O modelo de otimização robusta permite ao decisor incorporar a sua atitude em relação ao risco de forma muito simples. Assim, o decisor sabe exatamente qual é o pior caso teórico, diferente do modelo de programação estocástica, que deve ser resolvido a priori para se saber o pior caso.

3. iii

   O modelo de otimização robusta é mais tratável computacionalmente, pois os equivalentes robustos mantêm a complexidade computacional do problema determinístico, diferentemente dos modelos estocásticos. Porém, vale a ressalva de que conjuntos de incertezas não poliédricos (elipsoidais, por exemplo) podem tornar o equivalente robusto computacionalmente intratável.

4. iv

   O modelo de otimização robusta com poucos parâmetros incertos é, em geral, muito conservador, no sentido de que é preciso garantir quase o pior caso (modelo de Soyster) para assegurar uma pequena probabilidade de violação das restrições.

5. v

   Algumas análises próprias da otimização robusta podem não fazer sentido se a solução ótima não é determinada. Por exemplo, as probabilidades de violação das restrições referem-se a uma situação na qual a solução ótima do problema robusto existe e pode ser determinada. Caso contrário, os limitantes probabilísticos podem ser inócuos.

6. vi

   Não é simples determinar os equivalentes robustos quando várias fontes de incerteza afetam simultaneamente o problema. Nesse caso, é preciso interpretar as correlações entre os diversos parâmetros incertos, de modo que a versão robusta não seja muito conservadora.

7. vii

   A seleção dos budgets de incerteza no modelo de otimização robusta com demanda incerta não é trivial.

5 Sumário e considerações finais

Neste trabalho, foram comparados dois modelos de planejamento da produção em empresas moveleiras com incerteza na demanda, um modelo de programação estocástica de dois estágios com aversão ao risco e um modelo de otimização robusta com intervalo de incerteza poliédrico-intervalar. Foram analisadas três importantes características dos modelos: o aumento no valor ótimo da função objetivo, o nível de serviço e o tempo de execução computacional. Os resultados confirmaram a importância de escolher budgets de incerteza plausíveis, de modo a obter um tradeoff aceitável entre custo e robustez. Verificou-se também que o modelo equivalente robusto é mais fácil de ser resolvido do que a versão estocástica, o que é especialmente importante quando a versão determinística já apresenta dificuldade de resolução. Cabe ressaltar que a comparação feita neste trabalho é ilustrativa. Quais modelos ou técnicas são mais adequados em cada aplicação é uma questão mais difícil de ser respondida e depende do contexto da aplicação e do tomador de decisões. Mesmo a comparação apresentada aqui tem suas limitações, pois apenas um exemplar foi analisado, e considerando-se somente uma fonte de incerteza, no caso as demandas dos produtos. Para recomendações mais gerais sobre quais metodologias são mais indicadas no tratamento do problema combinado CLC sob incertezas, seria necessário realizar um planejamento de experimentos com diversos exemplares e com diferentes características, assim como considerar conjuntamente a presença de mais de uma fonte de incerteza (por exemplo, outros parâmetros além das demandas dos produtos). Convém salientar que mesmo que os resultados aqui apresentados sejam referentes a uma aplicação particular na indústria de móveis, é possível estender o paradigma apresentado a outros problemas correlatos.

Agradecimentos

O primeiro autor agradece as bolsas de pesquisa à FAPESP (processos 06/05054-8 e 2013/08303-2) e ao CNPq(306237/2014-8), o apoio financeiro.

• Suporte financeiro: FAPESP (processos 06/05054-8 e 2013/08303-2) e CNPq (306237/2014-8).

Referências

• Ahmed, S. (2006). Convexity and decomposition of mean-risk stochastic programs. Mathematical Programming Series A, 106(3), 433-446. http://dx.doi.org/10.1007/s10107-005-0638-8.
  » http://dx.doi.org/10.1007/s10107-005-0638-8
• Alem, D. J. (2011). Programação estocástica e otimização robusta no planejamento da produção de empresas moveleiras (Tese de doutorado). Universidade de São Paulo, São Carlos.
• Alem, D., & Morabito, R. (2012). Production planning in furniture settings via robust optimization. Computers & Operations Research, 39(2), 139-150. http://dx.doi.org/10.1016/j.cor.2011.02.022.
  » http://dx.doi.org/10.1016/j.cor.2011.02.022
• Alem, D., & Morabito, R. (2013a). O problema combinado de planejamento da produção e corte de estoque sob incertezas: aplicação em fábricas de móveis de pequeno porte. Gestão & Produção, 20(1), 111-133. http://dx.doi.org/10.1590/S0104-530X2013000100009.
  » http://dx.doi.org/10.1590/S0104-530X2013000100009
• Alem, D., & Morabito, R. (2013b). Risk-averse two-stage stochastic programs in furniture plants. OR-Spektrum, 35(4), 773-806. http://dx.doi.org/10.1007/s00291-012-0312-5.
  » http://dx.doi.org/10.1007/s00291-012-0312-5
• Alonso-Ayuso, A., Carvallo, F., Escudero, L. F., Guignard, M., Pi, J., Puranmalka, R., & Weintraub, A. (2014). Medium range optimization of copper extraction planning under uncertainty in future copper prices. European Journal of Operational Research, 233(3), 711-726. http://dx.doi.org/10.1016/j.ejor.2013.08.048.
  » http://dx.doi.org/10.1016/j.ejor.2013.08.048
• Amorim, P., Alem, D., & Almada-Lobo, B. (2013). Risk management in production planning of perishable goods. Industrial & Engineering Chemistry Research, 52(49), 17538-17553. http://dx.doi.org/10.1021/ie402514c.
  » http://dx.doi.org/10.1021/ie402514c
• Barreiros, D. A., Alem, D., & Toso, E. A. V. (2013). Planejamento agregado na indústria de nutrição animal sob incertezas. Produção. No prelo.
• Ben-Tal, A., Golany, B., Nemirovski, A., & Vial, J. (2005). Supplier-retailer flexible commitments contracts: a robust optimization approach. Manufacturing & Service Operations Management, 7(3), 248-271. http://dx.doi.org/10.1287/msom.1050.0081.
  » http://dx.doi.org/10.1287/msom.1050.0081
• Ben-Tal, A., & Nemirovski, A. (2000). Robust solutions of linear programming problems contaminated with uncertain data. Mathematical Programming Series B, 88(3), 411-424. http://dx.doi.org/10.1007/PL00011380.
  » http://dx.doi.org/10.1007/PL00011380
• Bertsimas, D., & Sim, M. (2003). Robust discrete optimization and network flows. Mathematical Programming, 98(1-3), 43-71. http://dx.doi.org/10.1007/s10107-003-0396-4.
  » http://dx.doi.org/10.1007/s10107-003-0396-4
• Bertsimas, D., & Sim, M. (2004). The price of robustness. Operations Research, 52(1), 35-53. http://dx.doi.org/10.1287/opre.1030.0065.
  » http://dx.doi.org/10.1287/opre.1030.0065
• Bertsimas, D., & Thiele, A. (2006). A robust optimization approach to inventory theory. Operations Research, 54(1), 150-168. http://dx.doi.org/10.1287/opre.1050.0238.
  » http://dx.doi.org/10.1287/opre.1050.0238
• Birge, J. R., & Louveaux, F. (1997). Introduction to stochastic programming. New York: Springer.
• Bohle, C., Maturana, S., & Vera, J. (2009). A robust optimization approach to wine grape harvesting scheduling. European Journal of Operational Research, 200(1), 245-252. http://dx.doi.org/10.1016/j.ejor.2008.12.003.
  » http://dx.doi.org/10.1016/j.ejor.2008.12.003
• Gramani, M., & França, P. (2006). The combined cutting stock and lot-sizing problem in industrial process. European Journal of Operational Research, 174(1), 509-521. http://dx.doi.org/10.1016/j.ejor.2004.12.019.
  » http://dx.doi.org/10.1016/j.ejor.2004.12.019
• ILOG. (2008). ILOG CPLEX 11.0: user's manual and reference manuals. Recuperado em junho de 2009, de http://www.ilog.com/products/cplex/
  » http://www.ilog.com/products/cplex/
• Kall, P., & Wallace, S. (1994). Stochastic programming. New York: Wiley.
• Li, Z., & Ierapetritou, M. G. (2008). Robust optimization for process scheduling under uncertainty. Industrial & Engineering Chemistry Research, 47(12), 4148-4157. http://dx.doi.org/10.1021/ie071431u.
  » http://dx.doi.org/10.1021/ie071431u
• Markowitz, H. (1959). Portfolio selection: efficient diversification of investments. New York: John Wiley & Sons.
• McCarl, B., Meeraus, A., Eijk, P., Bussieck, M., Dirkse, S., & Steacy, P. GAMS user guide. Recuperado em março de 2010, de http://www.gams.com/dd/docs/bigdocs/gams2002/mccarlgamsuserguide.pdf
  » http://www.gams.com/dd/docs/bigdocs/gams2002/mccarlgamsuserguide.pdf
• Mulvey, J., Vanderbei, R., & Zenios, S. (1995). Robust optimization of large-scale systems. Operations Research, 43(2), 264-281. http://dx.doi.org/10.1287/opre.43.2.264.
  » http://dx.doi.org/10.1287/opre.43.2.264
• Munhoz, J. R., & Morabito, R. (2013). Uma abordagem de otimização robusta no planejamento agregado de produção na indústria cítrica. Produção, 23(2), 422-435.
• Paiva, R. P. O., & Morabito, R. (2014). Planejamento hierárquico da produção em usinas de açúcar e álcool: modelagem de otimização robusta. Produção, 24(3),644-663.
• Palma, C. D., & Nelson, J. D. (2009). A robust optimization approach protected harvest scheduling decisions against uncertainty. Canadian Journal of Forest Research, 39(2), 342-355. http://dx.doi.org/10.1139/X08-175.
  » http://dx.doi.org/10.1139/X08-175
• Rangel, S., & Figueiredo, A. G. (2008). O problema de corte de estoque em indústrias de móveis de pequeno e médio portes. Pesquisa Operacional, 28(3), 451-472. http://dx.doi.org/10.1590/S0101-74382008000300004.
  » http://dx.doi.org/10.1590/S0101-74382008000300004
• Rockafellar, R. T., & Uryasev, S. (2000). Optimization of conditional value-at-risk. Journal of Risk, 2, 21-41.
• Shapiro, A., Dentcheva, D., & Ruszczynski, A. (2009). Lectures on stochastic programming: modeling and theory. Philadelphia: SIAM. http://dx.doi.org/10.1137/1.9780898718751.
  » http://dx.doi.org/10.1137/1.9780898718751
• Soyster, A. (1973). Convex programming with set-inclusive constraints and applications to inexact linear programming. Operations Research, 21(5), 1154-1157. http://dx.doi.org/10.1287/opre.21.5.1154.
  » http://dx.doi.org/10.1287/opre.21.5.1154
• Takriti, S., & Ahmed, S. (2004). On robust optimization of two-stage systems. Mathematical Programming Series A, 99(1), 109-126. http://dx.doi.org/10.1007/s10107-003-0373-y.
  » http://dx.doi.org/10.1007/s10107-003-0373-y
• Toso, E. A. V., & Alem, D. (2014). Effective location models for sorting recyclables in public management. European Journal of Operational Research, 234(3), 839-860. http://dx.doi.org/10.1016/j.ejor.2013.10.035.
  » http://dx.doi.org/10.1016/j.ejor.2013.10.035
• Wei, C., Li, Y., & Cai, X. (2011). Robust optimal policies of production and inventory with uncertain returns and demand. International Journal of Production Economics, 134(2), 357-367. http://dx.doi.org/10.1016/j.ijpe.2009.11.008.
  » http://dx.doi.org/10.1016/j.ijpe.2009.11.008
• Zanjani, M. K., Ait-Kadi, D., & Nourelfath, M. (2009). Robust production planning in a manufacturing environment with random yield: A case in sawmill production planning. European Journal of Operational Research, 201(3), 882-891. http://dx.doi.org/10.1016/j.ejor.2009.03.041.
  » http://dx.doi.org/10.1016/j.ejor.2009.03.041

Datas de Publicação

• Publicação nesta coleção
  29 Set 2015
• Data do Fascículo
  Jul-Sep 2015

Histórico

• Recebido
  22 Jan 2014
• Aceito
  12 Maio 2015