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