Lab2.Rmd

---
title: Лабораторная работа 2. Линейное программирование
author: "Автор"
date: "10/02/2024"
output:
  pdf_document:
    latex_engine: lualatex
    includes:
      in_header: preamble.tex
  html_document:
    df_print: paged
  word_document: default
editor_options: 
  markdown: 
    wrap: sentence
---

# Графическое решение

$$\begin{cases}
      \max{F(x)=2x_1+6x_2}\\
      x_1+4x_2\leq2048\\
      2x_1+x_2\leq2048\\
      x2\leq480\\
      x1,x2\geq0
\end{cases}\,$$

Задача 1 найти область удовлетворяющую ограничениям.
Сразу можно сказать, что рассматривается первая четверть из-за ограничения на положительность переменных.
Определим функцию для построения замкнутой области на графике используя библиотеку **plotly**.

```{r echo=FALSE}
library(plotly)

add_region <- function(fig, xy, name = "region", color) {
  add_trace(
    fig,
    x = xy[, 1],
    y = xy[, 2],
    type = "scatter",
    fill = "toself",
    fillcolor = color,
    hoveron = "points+fills",
    marker = list(color = "red"),
    line = list(color = color),
    text = apply(xy, 1, paste, collapse = ","),
    hoverinfo = "text",
    name = TeX(name)
  )
}
```

```{r echo=FALSE}
library(plotly)

plot_ly(mode = "lines+markers") |>
  add_region(
    xy = rbind(c(0, 0), c(0, 2048 / 4), c(2048, 0)),
    color = "rgba(255, 212, 96, 0.5)",
    name = "x_1+4x_2\\leq2048") |>
  add_region(
    xy = rbind(c(0, 0), c(0, 2048), c(2048 / 2, 0)),
    color = "rgba(168, 216, 234, 0.5)",
    name = "2x_1+x_2\\leq2048") |>
  add_region(
    xy = rbind(c(0, 0), c(0, 480), c(2048, 480), c(2048, 0)),
    color = "rgba(54, 216, 234, 0.5)",
    name = "x_2\\leq480") |>
  config(mathjax = "cdn")
```

Итоговая область поиска:

```{r echo=FALSE}
plot_ly(mode = "lines+markers") |>
  add_region(
    xy = rbind(c(0, 0), c(0, 480), c(128, 480), c(6144 / 7, 2048 / 7), c(2048 / 2, 0)),
    name = "region",
    color = "rgba(168, 216, 234, 0.5)") |>
  config(mathjax = "cdn")
```

Теперь можно добавить график уровня функции и совместить его с графиком ограничений

```{r echo=FALSE}
x <- seq(0, 1100, length.out = 100)
y <- seq(0, 500, length.out = 100)
z <- outer(x, y, FUN = \(x1, x2) {
  2 * x1 + 6 * x2
}) |> t()

plot_ly(x = x, y = y, z = z, type = "contour") |>
  add_region(
    xy = rbind(c(0, 0), c(0, 480), c(128, 480), c(6144 / 7, 2048 / 7), c(2048 / 2, 0)),
    color = "rgba(168, 216, 234, 0.5)") |>
  config(mathjax = "cdn")
```

Следовательно наилучшее значение функции в точке $(x=\frac{6144}{7},y=\frac{2048}{7})$

# Численное решение

Существует два интерфейса для работы с этой библиотекой.
Более высокоуровненый lpSolve и низкоуровневый lpSolveAPI.
Рассмотрим оба

## lpSolve

Требуется задать коэффициенты целевой функции.
$c(a_1,a_2,a_3,...,a_n)\implies F(\vec{x})=a_1x_1+a_2x_2+a_3x_3+...a_nx_n$.

```{r}
library(lpSolve)
Fun <- c(2, 6)
```

Ограничения задаются в виде матрицы коэффициентов

```{r}
A <- rbind(c(1, 4), c(2, 1), c(0, 1))
```

Правые части и знаки ограничений задаются отдельно

```{r}
B <- c(2048, 2048, 480)
CD <- c("<=", "<=", "<=")
```

Для решения вызываетя процедура **lp**

```{r}
optimum <- lp(
  direction = "max",
  objective.in = Fun,
  const.mat = A,
  const.dir = CD,
  const.rhs = B,
  compute.sens = TRUE)
optimum
```

Макисмальное значение целевой функции достигается при следующих значениях переменных

```{r}
optimum$solution
```

## lpSolveAPI

Создать объект модели задачи линейного программирования с 2 переменными и 0 ограничений (будут добавлены позднее).

```{r}
library(lpSolveAPI)
lprob <- make.lp(0, 2)
```

Далее выбирается тип оптимизации, в данном случае максимизация.
Вывод функции перенаправляется в переменную заглушку, чтобы отменить печать информации на экран.

```{r}
stub <- lp.control(lprob, sense = "max")
rm(stub)
```

Далее задаются ограничения.
Коэффициенты в функция ограничения задаются с помощью $c(a_1,a_2,a_3,...,a_n)\implies a_1x_1+a_2x_2+a_3x_3+...a_nx_n$.
По умолчанию на переменные накладывается условие $\forall{i}: x_i\geq0$

```{r}
add.constraint(lprob, c(1, 4), "<=", 2048)
add.constraint(lprob, c(2, 1), "<=", 2048)
add.constraint(lprob, c(0, 1), "<=", 480)
```

Задание целевой функции.
Объект модели можно конструировать в произвольном порядке

```{r}
set.objfn(lprob, c(2, 6))
lprob
```

Далее необходимо решить систему

```{r}
status_code <- solve(lprob)
```

Наилучшее значение целевой функции

```{r}
get.objective(lprob)
```

Значение переменных, при которых дистигнуто максимальное значение целевой функции

```{r}
get.variables(lprob)
```

## Анализ чувствительности

Анализ чувствительности применим только к задаче в действительных числах.
Исходная целевая функция $F=2x_1+6x_2$.
Для анализа чувствительности целевой функции используется функция **get.sensitivity.obj**.

```{r}
get.sensitivity.obj(lprob)
```

Полученные значения означают, что при измениении значения $x_1\in(1.5,4)$ или $x_2\in(4,32)$ точка соовтетствующая лучшему значению целевой функции не изменится.

Для анализа чувствительности ограничений применяется функция **get.sensitivity.rhs**.

```{r}
get.sensitivity.rhs(lprob)
```

Рассмотрим первое ограничение $x_1+4x_2\leq2048$.
Значение $\xi_1=1.4285$ показывает насколько изменится целевая функция, если значение ограничения увеличится на одну единицу.
Значение $xi$ ненулевое, только если ограничение активно.
$\xi$ можно интерпретировать как частную произовдную целевой фунцкии по правой части неравенстра ограничения.
$\xi_1=\frac{\partial F}{\partial B_1}$.
Следовательно при изменении ограничения $xi_i$ в пределах от $xi_from$ до $xi_{till}$ значение целевой функции изменится как $F=F_0+\Delta B_i \xi_i$

Значение $\xi_i$ дает очень хорошее представление о том, сколько стоит это ограничение.
Если значение $\xi_i$ очень велико, то это ограничение очень влияет на целевую функцию, и если вы сможете немного изменить его, то решение будет намного лучше.
Также имеет значение знак $\xi_i$.
Положительное значение означает, что по мере увеличения ограничения значение целевой функции будет больше, а по мере того, как оно становится более отрицательным, значение целевой функции будет меньше.

Аналогично для *lpSolve* используются следующие функции для анализа коэффициентов целевой функции

```{r}
optimum$sens.coef.from
optimum$sens.coef.to
```

и для анализа ограничений

```{r}

```

# Особые случаи применеия симплекс метода

## Альтернативные решения

$$\begin{cases}
      \max{F(x)=x_1+2x_2}\\
      x_1\leq80\\
      5x_1+6x_2\leq600\\
      x_1+2x_2\leq160\\
      x1,x2\geq0
\end{cases}\,$$

Изобразим область поиска соответствующую ограничениям

```{r echo=FALSE}
plot_ly(mode = "lines+markers") |>
  add_region(
    xy = rbind(c(0, 100), c(120, 0), c(0, 0)),
    name = "5x_1+6x_2\\leq600",
    color = "rgba(255, 216, 234, 0.5)") |>
  add_region(
    xy = rbind(c(0, 80), c(160, 0), c(0, 0)),
    name = "x_1+2x_2\\leq160",
    color = "rgba(168, 216, 234, 0.5)") |>
  add_region(
    xy = rbind(c(0, 0), c(0, 100), c(80, 100), c(80, 0)),
    name = "x_1\\leq80",
    color = "rgba(168, 216, 234, 0.5)") |>
  config(mathjax = "cdn")
```

Теперь можно добавить график уровня функции и совместить его с графиком ограничений

```{r echo=FALSE}
x <- seq(0, 80, length.out = 100)
y <- seq(0, 80, length.out = 100)
z <- outer(x, y, FUN = \(x1, x2) {
  x1 + 2 * x2
}) |> t()

plot_ly(x = x, y = y, z = z, type = "contour") |>
  add_region(
    xy = rbind(c(0, 0), c(0, 80), c(60, 50), c(80, 100 / 3), c(80, 0)),
    name = "region",
    color = "rgba(168, 216, 234, 0.5)") |>
  config(mathjax = "cdn")
```

Из графика следует то, что прямая постоянного уровня целевой функции $F(x)=x_1+2x_2$ совпала с прямой задающей ограничение $x_1+2x_2\leq160$.
Значит все точки принадлежащие прямой проходящей через координаты $P(0,80)$ и $Q(60,50)$ будут соответствовать оптимальному решению.

```{r}
library(lpSolve)
Fun <- c(1, 2)
A <- rbind(c(0, 1), c(5, 6), c(1, 2))
B <- c(80, 600, 160)
CD <- c("<=", "<=", "<=")

optimum <- lp(
  direction = "max",
  objective.in = Fun,
  const.mat = A,
  const.dir = CD,
  const.rhs = B,
  compute.sens = TRUE)

optimum
optimum$solution
```

Процедура поиска решения находит только одну точку.

```{r}
optimum$sens.coef.from
optimum$sens.coef.to
```

Анализ чувствительности показывает, что активно только ограничение $x_1+2x_2\leq160$, а соответствующее значенеи $\xi=1$

```{r}
optimum$duals
optimum$duals.from
optimum$duals.to
```

## Отсутствие допустимых решений

В данном случае условия ограничений не совместимы.

$$\begin{cases}
      \max{F(x)=200x_1+300x_2}\\
      2x_1+3x_2\geq1200\\
      x_1+x_2\leq400\\
      2x_1+1.5x_2\geq900\\
      x1,x2\geq0
\end{cases}\,$$

```{r echo=FALSE}
plot_ly(mode = "lines+markers") |>
  add_region(
    xy = rbind(c(600, 0), c(0, 400), c(0, 600), c(600, 600)),
    name = "2x_1+3x_2\\geq1200",
    color = "rgba(255, 212, 255, 0.5)") |>
  add_region(
    xy = rbind(c(0, 0), c(400, 0), c(0, 400)),
    name = "x_1+x_2\\leq400",
    color = "rgba(255, 212, 96, 0.5)") |>
  add_region(
    xy = rbind(c(0, 600), c(450, 0), c(600, 0), c(600, 600)),
    name = "2x_1+1.5x_2\\geq900",
    color = "rgba(255, 255, 96, 0.5)") |>
  config(mathjax = "cdn")
```

```{r}
library(lpSolve)
Fun <- c(200, 300)
A <- rbind(c(2, 3), c(1, 2), c(2, 1.5))
B <- c(1200, 400, 900)
CD <- c(">=", "<=", ">=")

optimum <- lp(
  direction = "max",
  objective.in = Fun,
  const.mat = A,
  const.dir = CD,
  const.rhs = B,
  compute.sens = TRUE)

optimum
optimum$solution
```

Процедура поиска решения завершается с ошибкой

## Неограниченность решений

It is a solution whose objective function is infinite.
If the feasible region is unbounded then one or more decision variables will increase indefinitely without violating feasibility, and the value of the objective function can be made arbitrarily large.
Consider the following model:

$$\begin{cases}
      \max{F(x)=40x_1+60x_2}\\
      2x_1+x_2\geq70\\
      x_1+x_2\geq40\\
      x_1+3x_2\geq90\\
      x1,x2\geq0
\end{cases}\,$$

```{r echo=FALSE}
plot_ly(mode = "lines+markers") |>
  add_region(
    xy = rbind(c(0, 40), c(0, 35), c(70, 0), c(90, 0), c(90, 40)),
    name = "2x_1+x_2\\geq70",
    color = "rgba(255, 212, 255, 0.3)") |>
  add_region(
    xy = rbind(c(0, 40), c(40, 0), c(90, 0), c(90, 40)),
    name = "x_1+x_2\\leq40",
    color = "rgba(255, 212, 96, 0.3)") |>
  add_region(
    xy = rbind(c(0, 40), c(0, 30), c(90, 0), c(90, 40)),
    name = "x_1+3x_2\\geq90",
    color = "rgba(255, 255, 96, 0.3)") |>
  config(mathjax = "cdn")
```

Объединяя добавляя контурный график

```{r echo=FALSE}
x <- seq(0, 90, length.out = 100)
y <- seq(0, 40, length.out = 100)
z <- outer(x, y, FUN = \(x1, x2) {
  40 * x1 + 60 * x2
}) |> t()

plot_ly(x = x, y = y, z = z, type = "contour") |>
  add_region(
    xy = rbind(c(15, 40), c(24, 22), c(90, 0), c(90, 40)),
    name = "region",
    color = "rgba(168, 216, 234, 0.5)") |>
  config(mathjax = "cdn")
```

Ограничения заданные в этом примере не ограничивают область поиска в направлениии $x\to\infty$ и $y\to\infty$.
Целевая функция неограниченно возрастает в этом направлении, следовательно и решение неограниченно.

```{r}
library(lpSolve)
Fun <- c(40, 60)
A <- rbind(c(2, 1), c(1, 1), c(1, 3))
B <- c(70, 40, 90)
CD <- c(">=", ">=", ">=")

optimum <- lp(
  direction = "max",
  objective.in = Fun,
  const.mat = A,
  const.dir = CD,
  const.rhs = B,
  compute.sens = TRUE)

optimum
optimum$solution
```

Процедура поиска решения завершается с ошибкой 3: *UNBOUNDED (3) The model is unbounded*.