2  Графичен метод

Покажи 
# %pip install plotly --quiet

from plotly import graph_objects as go

2.1 Лекция

2.2 Допълнителни материали

2.3 Линейно програмиране в Excel

2.4 Задача: планиране на производство

Бутиково кафене предлага два вида еспресо: Супер еспресо и Делукс еспресо. За приготвянето на един килограм от първия вид еспресо са необходими по равни части арабика и робуста, а рецептата за Делукс предвижда смес от арабика и робуста в пропорция 1 към 3. Доставчиците са готови да осигурят 120 кг арабика и 160 кг. робуста. Заведението знае, че няма да може да продаде повече от 150 кг. Делукс еспресо. От всеки продаден килограм Супер еспресо заведението печели 40 EUR, докато печалбата от килограм Делукс възлиза на 50 EUR.

Колко от двата типа еспресо ще препоръчате на кафенето да произведе, за да максимизира печалбата си?

Целеви променливи:

\begin{align*} & x: \text{ Супер еспресо (кг.)}\\ & y: \text{ Делукс еспресо (кг.)} \end{align*}

\max z = 40 x + 50 y \text{ (целева функция)}

\begin{align*} 0.5 x + 0.25 y & \leq 120 \text{ (арабика)}\\ 0.5 x + 0.75 y & \leq 160 \text{ (робуста)} \\ 0 \cdot x + y & \leq 150 \text{ (търсене Делукс)}\\ x & \geq 0 \\ y & \geq 0 \end{align*}

2.5 Допустимо множество

Както и в предишната задача ще изобразим графично допустимото множество, като начертаем правите, към всяко от петте неравенства:

\begin{align} 0.5 x & + 0.25 y & = & 120 \\ 0.5 x & + 0.75 y & = & 160 \\ 0 \cdot x & + y & = & 150 \end{align}

Първо ще пресметнем пресечните точки на трите прави с осите x и y?

  • Права арабика: (0, 120 / 0.25 = 480), (120 / 0.5 = 240, 0)
  • Права робуста: (0, 160 / 0.75), (160 / 0.5, 0)
  • Права търсене на Делукс: (0, 150), (100, 150) Тази права е успоредна на оста x.
Покажи 
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display

fig, ax = plt.subplots()

ax.plot([0, 240], [480, 0], label=r"$0.5 x + 0.25 y = 120$ (арабика)")
ax.plot([0, 320], [213.33, 0], label=r"$0.5 x + 0.75 y = 160$ (робуста)")
ax.plot([0, 400], [150, 150], label=r"$y = 150$ (търсене на Делукс)")

ax.set_xlim((0, 400))
ax.set_ylim((0, 500))
ax.set_xlabel(r'$x$ (Супер, кг.)')
ax.set_ylabel(r'$y$ (Делукс, кг.)')
ax.legend(loc=0)
Фигура 2.1: Прави на трите ограничения
Покажи 
vertices = [
    [0, 0],
    [0, 150],
    [95, 150],
    [200, 80],    
    [240, 0]
]

vertices_x = [v[0] for v in vertices]
vertices_y = [v[1] for v in vertices]

for v in vertices:
    ax.annotate(
        f"({v[0]}, {v[1]})",
        (v[0], v[1]),
        textcoords="offset points",
        xytext=(0, 5)
    )

ax.fill(vertices_x, vertices_y, color='grey', alpha=0.3)

display(fig)
Фигура 2.2: Допустимо множество

2.6 Целева функция и оптимален план

За да определим оптималния план графично ще начертаем прави, съответстващи на различни нива на печалба.

Целевата функция зависи от две променливи, което означава, че за да я изобразим графично ще трябва да използваме трето измерение. Фигура 2.3 показва графиката на целевата функция в тримерно пространство, както и нейната проекция в двумерно пространство.

Покажи 
# Create a 3d scatter plot for the vertices of the feasible region with plotly

fig = go.Figure(data=[go.Scatter3d(
    x=vertices_x,
    y=vertices_y,
    z=[0] * len(vertices_x),
    mode='markers+text',
    marker=dict(size=2, color='red'),
    text=[f"({v[0]}, {v[1]}, 0)" for v in vertices],
    textposition="top center"
)])

# Draw lines connecting the vertices

for i in range(len(vertices)):
    x_line = [vertices[i][0], vertices[(i + 1) % len(vertices)][0]]
    y_line = [vertices[i][1], vertices[(i + 1) % len(vertices)][1]]
    z_line = [0, 0]
    fig.add_trace(go.Scatter3d(
        x=x_line,
        y=y_line,
        z=z_line,
        mode='lines',
        line=dict(color='blue', width=2)
    ))

fig.update_layout(
    scene=dict(
        xaxis_title=r"Super (kg)",
        yaxis_title=r"Deluxe (kg)",
        zaxis_title=r"Profit (Euro)",
        zaxis=dict(range=[0, 12000]),
        aspectmode='cube'
    ),
    showlegend=False
)

# Draw the plane of the objective function

X1 = np.linspace(0, 250, 20)
X2 = np.linspace(0, 250, 20)

X1, X2 = np.meshgrid(X1, X2)
Z = 40 * X1 + 50 * X2

fig.add_trace(go.Surface(
    x=X1,
    y=X2,
    z=Z,
    opacity=0.5,
    colorscale='Viridis'
))

# Draw a plane at level z = 6000

z_plane_1 = 6000

fig.add_trace(go.Surface(
    x=X1,
    y=X2,
    z=np.full(X1.shape, z_plane_1),
    opacity=0.2,
    colorscale='Viridis',
    showscale=False
))

# Draw a line at the intersection of the objective plane and the z = 10000 plane
intersection_x = np.linspace(0, 250, 20)
intersection_y_1 = (z_plane_1 - 40 * intersection_x) / 50

fig.add_trace(go.Scatter3d(
    x=intersection_x,
    y=intersection_y_1,
    z=[z_plane_1] * len(intersection_x),
    mode='lines',
    opacity=0.4,
    line=dict(color='red', width=3),
    name='Intersection'
))

# Draw the intersection line on the xy plane

fig.add_trace(go.Scatter3d(
    x=intersection_x,
    y=intersection_y_1,
    z=[0] * len(intersection_x),
    mode='lines',
    opacity=0.4,
    line=dict(color='red', width=2),
    name='Intersection'
))
fig.update_layout(
    title="Profit Function, Feasible Region and Isolines at 6000 and 12000 Euro",
    scene=dict(
        xaxis_title=r"Super (kg)",
        yaxis_title=r"Deluxe (kg)",
        zaxis_title=r"Profit (Euro)",
        zaxis=dict(range=[0, 13000]),
        yaxis=dict(range=[0, 250]),
        xaxis=dict(range=[0, 250]),
        aspectmode='cube'
    ),
    showlegend=False
)

# Draw a plane at the optimal solution z = 12000

z_plane_2 = 12000

fig.add_trace(go.Surface(
    x=X1,
    y=X2,
    z=np.full(X1.shape, z_plane_2),
    opacity=0.2,
    colorscale='Viridis',
    showscale=False
))

intersection_x = np.linspace(0, 250, 20)
intersection_y = (z_plane_2 - 40 * intersection_x) / 50

# Draw the intersection line on the xy plane

fig.add_trace(go.Scatter3d(
    x=intersection_x,
    y=intersection_y,
    z=[z_plane_2] * len(intersection_x),
    mode='lines',
    opacity=0.4,
    line=dict(color='red', width=3),
    name='Intersection'
))
fig.add_trace(go.Scatter3d(
    x=intersection_x,
    y=intersection_y,
    z=[0] * len(intersection_x),
    mode='lines',
    opacity=0.4,
    line=dict(color='red', width=3),
    name='Intersection'
))

fig.show()
Фигура 2.3: The profit function, the feasible region and the isolines at 6000 and 12000 Euro.

Нека да начертаем тези прави на графиката на допустимото множество. Фигура 2.4 показва графиката на допустимото множество и правите на печалбата от 3000 EUR, 50000 EUR и 12000 EUR

Покажи 
fig, ax = plt.subplots()

ax.fill(vertices_x, vertices_y, color='grey', alpha=0.3)
for v in vertices:
    ax.annotate(
        f"({v[0]}, {v[1]})",
        (v[0], v[1]),
        textcoords="offset points",
        xytext=(0, 5)
    )

ax.plot([0, 3000 / 40], [3000 / 50, 0], label=r'$40 x + 50 y = 3000$')
ax.plot([0, 5000 / 40], [5000 / 50, 0], label=r'$40 x + 50 y = 5000$')
ax.plot([0, 12000 / 40], [12000 / 50, 0], label=r'$40 x + 50 y = 12000$')

ax.quiver(0, 0, 40, 50, angles='xy', scale_units='xy', scale=1, color='grey', label='Нормален вектор')

ax.set_xlim((0, 260))
ax.set_ylim((0, 260))

ax.legend(loc=0)
plt.show()
Фигура 2.4: Оптимален план

Във Фигура 2.4 е показан векторът (40, 50)^T, който е перпендикулярен на всички прави на печалбата (които са успоредни една спрямо друга). Когато задачата се решава графично на ръка, този вектор е удобен, защото показва посоката на увеличение на печалбата (в нашия случай тя расте по-бързо в посока y, отколкото в посока x, защото коефициентът на y (50) е по-голям от коефициента на x (40)).

Защо този вектор е перпендикулярен на правите на печалба? Един лесен начин да видим това е да изразим правата в параметричен вид. За по-лесно нека фиксираме някакво ниво на печалба, например 5000 EUR. Всички комбинации от супер и делукс, които постигат точно това ниво на печалба лежат на права, дадена от следното уравнение:

40x + 50y = 5000

Това е т.н. координатно представяне на правата. Параметричното представяне използва една (произволна) начална точка и вектор, който показва посоката й. За да намерим началната точка, ще намерим пресечните точки на правата с двете оси. Тази права се пресича с оста x при x = 125 и с оста y при y = 100, т.е. две точки на правата са (125, 0) и (0, 100).

От тези две точки можем да изведем параметричното представяне на правата (за t \in \mathbb{R}) като вземем разликата между двете точки (вектори).

\begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} 125 \\ 0 \end{pmatrix} + t \left[ \begin{pmatrix} 0 \\ 100 \end{pmatrix}- \begin{pmatrix}125 \\ 0 \end{pmatrix} \right] = \begin{pmatrix} 125 \\ 0 \end{pmatrix} + t \begin{pmatrix}-125 \\ 100 \end{pmatrix}

Два вектора са перпердикулярни (ортогонални), ако скаларното им произведение е равно на 0. Скаларното произведение на два вектора е равно на сумата от произведенията на съответните им координати. Така получаваме

\begin{pmatrix} -125 & 100 \end{pmatrix} \begin{pmatrix} 40 \\ 50 \end{pmatrix} = (-125) \cdot 40 + 100 \cdot 50 = -5000 + 5000 = 0

Упражнение 2.1 (Oще една задача от същия вид) Предприятие произвежда два вида боя: екстериорна (за външни повърхности) и интериорна (за вътрешни помещения). Производството двете бои изисква три суровини: разтворител, багрило и смола. За производство на един тон екстериорна боя са нужни 6 тона разтворител, един тон багрило и 2 тона смола. За производство на един тон интериорна боя са нужни 4 тона разтворител, 2 тона багрило и 7 тона смола. На разположение са следните количества суровини: 24 тона разтворител, 6 тона багрило и 4 тона смола. Колко тона от всяка боя да произведе предприятието, за да постигне максимална печалба? От продажбата на един тон екстериорна боя предприятието печели 5 000 EUR, а от продажбата на един тон интериорна боя - 4 000 EUR

Таблица 2.1: Консумация на материали за производството на интериорна и екстериорна боя и печалба за продаден тон
Екстериорна боя Интериорна боя Дневна наличност на ресурс
Разтворител 6 4 24
Багрило 1 2 6
Смола 2 7 4
Печалба (1 000 EUR/тон) 5 4

Съставете математически модел, който да намери оптималния план за производство на боите.

Целеви променливи: \begin{align*} & x: \text{ Екстериорна боя (тонове)}\\ & y: \text{ Интериорна боя (тонове)} \end{align*}

\max z = 5 x + 4 y \text{ (целева функция)}

Ограничения:

\begin{align*} 6x + 4y & \leq 24 \text{ (разтворител)}\\ x + 2y & \leq 6 \text{ (багрило)}\\ 2x + 7y & \leq 4 \text{ (смола)}\\ x & \geq 0 \\ y & \geq 0 \end{align*}

Покажи 
import gurobipy as gp
from gurobipy import GRB

# Create a new model
m = gp.Model("Paint")
m.Params.LogToConsole = 0

# Create variables
p_exterior = m.addVar(lb=0, name="exterior")
p_interior = m.addVar(lb=0, name="interior")

m.setObjective(5 * p_exterior + 4 * p_interior, GRB.MAXIMIZE)

m.addConstr(6 * p_exterior + 4 * p_interior <= 24, "Solvent")
m.addConstr(1 * p_exterior + 2 * p_interior <= 6, "Dye")
m.addConstr(2 * p_exterior + 7 * p_interior <= 4, "Resin")

m.optimize()

print(f"Optimal solution: p_exterior = {p_exterior.x}, p_interior = {p_interior.x}")
Restricted license - for non-production use only - expires 2027-11-29
Set parameter LogToConsole to value 0
Optimal solution: p_exterior = 2.0, p_interior = 0.0