Como encontrar os limitantes superior e inferior de um direito contigente.

Question

Considere o seguinte mercado:

\(S = 4; J = 2, p = (1,1), X = \left[\begin{array}{cc}
0.8 & 1 \\
0.9 & 1 \\
1.2 & 1 \\
1.3 & 1
\end{array}\right]\)

Considere o problema de encontrar os limitantes inferior e superior para o direito contigente dado por \(\left[\begin{array}{cc}
-2&-1&3&5\end{array}\right]\):
a) Encontre a solução numérica usando a ferramenta de programação linear do python.
b) Apresente a solução gráfica do problema para o caso do limitante superior.

Answer 1

a)
O limitante superior é dado por \(q_u(z) = min \{p.h: Xh \ge z\}\) e o limitante inferior por \(q_l(z) = max \{p.h: Xh \le z\}\). Os valores serão encontrados utilizando a biblioteca scipy do python.

import numpy as np
from scipy.optimize import linprog
#definindo as variáveis do problema
p = np.array([1.0, 1.0])

X = np.array([[0.8, 1.0], [0.9,1.0], [1.2,1.0], [1.3, 1.0]])

z = np.array([[-2], [-1], [3], [5]])

A função de programação linear do scipy resolve o problema de programação linear da seguinte forma: \(min\{ph: Xh \le z\}\). Então, para resolver \(q_u\) multiplicaremos as restrições por \(-1\) e para resolver \(q_l\) multiplicaremos a função objetivo por \(-1\). Assim os problemas ficam de acordo com as definições da função.

Resolvendo \(q_u\):

res = linprog(p, A_ub = np.multiply(-1.0,X), b_ub = np.multiply(-1.0,z), bounds = ((None,None), (None,None)),options={"disp": True})
print(res)

Optimization terminated successfully.
Current function value: 0.800000
Iterations: 3
fun: 0.7999999999999985
message: 'Optimization terminated successfully.'
nit: 3
slack: array([0. , 0.4, 0.6, 0. ])
status: 0
success: True
x: array([ 14. , -13.2])

Resolvendo \(q_l\):

res2 = linprog(np.multiply(-1.0,p), A_ub = X, b_ub = z, bounds = ((None,None), (None,None)),method = 'simplex',options = {'disp': True})
print(res2)

Optimization terminated successfully.
Current function value: -0.333333
Iterations: 4
fun: -0.3333333333333331
message: 'Optimization terminated successfully.'
nit: 4
slack: array([0.33333333, 0.66666667, 0. , 0. ])
status: 0
success: True
x: array([ 13.33333333, -13. ])

Segue solução alternativa para \(q_l\) utilizando a biblioteca Pulp do python que permite construir o problema de programação linear sem forma definida:

from pulp import *

prob = LpProblem("limitante inferior", LpMaximize)

LpVariable("h", None, None)

h1 = LpVariable("h1", None, None)
h2 = LpVariable("h2", None,None)
#você "adiciona" as equações ao problema
prob += 1.0*h1 + 1.0*h2
prob += 0.8*h1 + 1.0*h2 <= -2.0
prob += 0.9*h1 + 1.0*h2 <= -1.0
prob += 1.2*h1 + 1.0*h2 <= 3.0
prob += 1.3*h1 + 1.0*h2 <= 5.0

prob.solve()

print("Status:", LpStatus[prob.status])

for v in prob.variables():
    print(v.name, "=", v.varValue)

Status: Optimal
h1 = 13.333333
h2 = -13.0

b)

Para a solução gráfica do problema, criei um gráfico interativo utilizando a biblioteca bokeh do python. Publiquei o gráfico na seguinte página. Na página existem instruções para utilizar. solução gráfica
(O gráfico está com um problema que o ponto de mínimo (0.8) do problema nem sempre aparece, se você ficar mexendo ali na vizinhança ele aparece. Quando eu identificar a causa, edito a publicação.)

Segue o código para criação do gráfico:
(Detalhe: Como o código exporta o gráfico para HTML, as ações que alteram o gráfico devem ser escritas em javascript)

import numpy as np
from bokeh.plotting import figure, output_file, show, output_notebook
from ipywidgets import interact
from bokeh.io import push_notebook, show, output_notebook
from bokeh.layouts import widgetbox,column
from bokeh.models.widgets import Slider
from bokeh.models import CustomJS, ColumnDataSource, Title

def eq1(x):
    return -2 - 0.8*x

def eq2(x):
    return -1 - 0.9*x

def eq3(x):
    return 3 - 1.2*x

def eq4(x):
    return 5 - 1.3*x

def eq_obj(x, c):
    return c - x

x = np.linspace(-28.0,28.0,1400)
x[0] = 14.0
y = eq_obj(x, 0)
source = ColumnDataSource(data=dict(x=x, y=y))
x_fac = []
y_fac = []

source2 = ColumnDataSource(data=dict(x_fac=x_fac, y_fac = y_fac))

output_file("grafico.html")

p = figure(tools = "pan,box_zoom,reset,save", y_range = [-15,15], x_range = [-15,15], title = "limitante superior", x_axis_label = "h1", y_axis_label = "h2",width=600,height=600, sizing_mode='scale_width')

p.line(x, eq1(x), legend = "restr1", line_color = "green")
p.line(x, eq2(x), legend = "restr2", line_color = "red")
p.line(x, eq3(x), legend = "restr3", line_color = "blue")
p.line(x, eq4(x), legend = "restr4", line_color = "orange")
p.line('x', 'y', source = source,  legend = "eq_obj", line_color = "purple", line_dash = "4 4")
p.line(np.linspace(-100,100,100),0, line_color = "black")
p.line(0,np.linspace(-100,100,100), line_color = "black")
p.circle('x_fac', 'y_fac', source = source2, legend = "pontos factíveis", color = "navy", alpha = 0.5)

slider = Slider(start = -20, end = 20, value = 0, step = .01, title = "curva de nível")


callback = CustomJS(args= dict(source=source, source2 = source2), code="""

    function eq1(x) {
        return -2 - 0.8*x 
    }
    function eq2(x) {
        return -1 - 0.9*x 
    }
    function eq3(x) {
        return 3 - 1.2*x 
    }
    function eq4(x) {
        return 5 - 1.3*x 
    }
    function eq_obj(x, c) {
        return c - x 
    }

    var data = source.data
    var data2 = source2.data

    //valor do slider
    var c = cb_obj.value;

    var x = data['x']
    var y = data['y']

    data2['x_fac'] = []
    data2['y_fac'] = [] 

    for (var i = 0; i < x.length; i++) {
        y[i] = eq_obj(x[i], c)
        if (eq1(x[i]) <= y[i] && eq2(x[i]) <= y[i] && eq3(x[i]) <= y[i] && eq4(x[i]) <= y[i] ) {
            data2['x_fac'].push(x[i])
            data2['y_fac'].push(y[i])
        }
    }
    source.change.emit()
    source2.change.emit()

""")

slider.js_on_change('value', callback)


show(column(p, widgetbox(slider)))

Para mais informações: bokeh

Comments

Ótima resposta, Marcio. Um comentário sobre o gráfico: ao dar zoom, vemos que os pontos da região factível se afastam um dos outros, isso porque o computador trabalha, por natureza, com valores discretos. É importante notar que, para o limitante superior, a região factível é toda área acima da reta amarela, até a interseção com a reta verde, e a partir daí, acima da reta verde.