Пусть игрок А обладает m возможными стратегиями А₁,…,А_m, а природа П может находиться в одном из n своих состояний П₁,…,П_n. Предполагается обычно, что игрок А в состоянии оценить результаты выбора им каждой из своих стратегий А_i, i=1,…,m, при каждом состоянии природы П_j, j=1,…,n, количественно выражающиеся действительными числами а_ij.

Эти числа, называемые выигрышами игрока А, можно записать в виде матрицы:

library(knitr)
mtx<-matrix(c(2,1,2,
            3,5,6,
            1,4,2,
            4,3,1.5),ncol=4)
rownames(mtx)<-paste0("A~",1:nrow(mtx),"~")
colnames(mtx)<-paste0("П~",1:ncol(mtx),"~")
kable(mtx)

(w<-c(0.5,0.1,0.2,0.2))

## [1] 0.5 0.1 0.2 0.2

sum(w)

## [1] 1

B<-mtx %*% w
colnames(B)<-"П~B~"
kable(B)

max(B)

## [1] 2.4

(w<-c(0.25,0.25,0.25,0.25))

## [1] 0.25 0.25 0.25 0.25

BL<-mtx %*% w
colnames(BL)<-"П~B~"
kable(BL)

max(BL)

## [1] 3.25

V<-as.matrix(apply(mtx,1,min))
colnames(V)<-"П~W~"
kable(V)

max(V)

## [1] 1.5

kable(HL<-cbind(B,V))

(w<-c(0.75,0.25))

## [1] 0.75 0.25

kable(HLM<-HL %*% w,col.names = "КХЛ")

max(HLM)

## [1] 2.1

M<-as.matrix(apply(mtx,1,max))
colnames(M)<-"П~M~"
kable(M)

max(M)

## [1] 6

kable(HV<-cbind(V,M))

(w<-c(0.5,0.5))

## [1] 0.5 0.5

kable(HVM<-HV %*% w,col.names = "КГ")

max(HVM)

## [1] 3.75

mtx1<- matrix(rep(apply(mtx,2,max),nrow(mtx)),byrow=T,ncol=ncol(mtx))-mtx
kable(mtx1)

S<-as.matrix(apply(mtx1,1,max))
colnames(S)<-"П~S~"
kable(S)

min(S)

## [1] 1

1. Скачиваем и устанавливаем на флешку платформу PortableApp.
2. Скачиваем и устанавливаем R Portable.
3. Скачиваем и устанавливаем RStudio Portable.

Все необходимое установлено. Теперь настраиваем. Для этого делаем следующие действия.

1. Запускаем приложение Start  в корне флеш-накопителя. Видим в трее новый значок, если на него нажать открывается окно.

2. Запускаем RStudioPortable. Видим вот такое окно:

Необходимо указать путь к R. Для этого нажимаем кнопку Browse и на флеш-накопителе выбираем папку F:\PortableApps\R-Portable\App\R-Portable\bin (в моём случае это диск F, у вас может быть другой). Приобретает вот такой вид (поскольку у меня 64-битная ОС).

Выбираем ту, что больше нравится, я выбрал 64-битную и нажимаем OK. Видим вот такое окно.

Соглашаемся и получаем запущенную RStudio

Теперь необходимо установить пакеты расширения, которые нам пригодятся в работе. Для этого полностью скопируем команду

source(url("https://tushavin.ru/RStudio/installall.R"))

и вставим после символа приглашения к вводу (значок «>» в левом окне), нажимаем клавишу Enter и идем пить чай, ну, или делаем что-либо другое на компьютере, поскольку все прекрасно работает и в фоновом режиме. Некоторое время будут писаться разные красные и черненькие буковки, не обращаем внимание, значит функционирует. После того как все закончится снова появится символ приглашения ввода.

3. Переходим к проверке работы. Открываем материал по линейному программированию. Внизу видим раздел «Решение с пакетом lpSolve» и прямо под ним код. Копируем его полностью, вставляем после приглашения, нажимаем Enter. Вставляем еще две команды:

result
result$solution

Результат должен совпасть с описанным в заметке.

Задача. Определить оптимальный маршрут из пункта 1 в пункт 10 по схеме маршрута движения.

# Здесь и далее граф строится с помощью команды dot из пакета graphviz

digraph ex01 {
                rankdir=LR;
                size="8,5"
                node [shape = box];
                "1" -> "2" [ label = "3" ];
                "1" -> "3" [ label = "7" ];
                "1" -> "4" [ label = "2" ];
                "2" -> "5" [ label = "9" ];
                "2" -> "6" [ label = "11" ];
                "3" -> "5" [ label = "5" ];
                "3" -> "6" [ label = "10" ];
                "3" -> "7" [ label = "7" ];
                "4" -> "6" [ label = "15" ];
                "4" -> "7" [ label = "13" ];
                "5" -> "8" [ label = "7" ];
                "5" -> "9" [ label = "5" ];
                "6" -> "8" [ label = "3" ];
                "6" -> "9" [ label = "4" ];
                "7" -> "8" [ label = "7" ];
                "7" -> "9" [ label = "1" ];
                "8" -> "10" [ label = "1" ];
                "9" -> "10" [ label = "4" ];
}

Каждый квадрат на схеме изображает один из населенных пунктов, которые для удобства пронумерованы.

Стоимость проезда из пункта i в пункт j обозначим через \(с_{ij}\) (обозначено на стрелках). Требуется определить такой путь из пункта 1 в пункт 10, общая стоимость которого является минимальной.

Решение без применения компьютера

Решение: Воспользуемся формулой реккурентных соотношений Беллмана $$f_n(i)=\min_j\{c_{ij}+f_{n-1}(j)\},\ n=\bar{{1,N}},$$

где N — количество этапов в решении; \(f_n(i)\) — стоимость, отвечающая стратегии минимальных затрат для пути от пункта i, если до конечного пункта остается n шагов; \(P_n(i)\) — решение, позволяющее достичь  \(f_n(i)\). Начинаем поиск от конечного пункта.

n=1

$$f_1(8)=c_{8,10}=1, P_1(8)=10;$$

$$f_1(9)=c_{9,10}=4, P_1(9)=10;$$

n=2

$$f_2(5)=\min\{c_{5,8}+f_1(8);c_{5,9}+f_1(9)\}=8, P_2(5)=8;$$

$$f_2(6)=\min\{c_{6,8}+f_1(8);c_{6,9}+f_1(9)\}=4, P_2(6)=8;$$

$$f_2(7)=\min\{c_{7,8}+f_1(8);c_{7,9}+f_1(9)\}=5, P_2(7)=9;$$

digraph ex02 {
                rankdir=LR;
                size="8,5"
                node [shape = box];
                "1" -> "2" [ label = "3" ];
                "1" -> "3" [ label = "7" ];
                "1" -> "4" [ label = "2" ];
                "2" -> "5" [ label = "9" ];
                "2" -> "6" [ label = "11" ];
                "3" -> "5" [ label = "5" ];
                "3" -> "6" [ label = "10" ];
                "3" -> "7" [ label = "7" ];
                "4" -> "6" [ label = "15" ];
                "4" -> "7" [ label = "13" ];            
# Упростили    
                "5" -> "8" [ label = "7" ];
                "6" -> "8" [ label = "3" ];
                "7" -> "9" [ label = "1" ];
                "8" -> "10" [ label = "1" ];
                "9" -> "10" [ label = "4" ];
}

n=3

$$f_3(2)=\min\{c_{2,5}+f_2(5);c_{2,6}+f_2(6)\}=15, P_3(2)=6;$$

$$f_3(3)=\min\{c_{3,5}+f_2(5);c_{3,6}+f_2(6);c_{3,7}+f_2(7)\}=12, P_3(3)=7;$$

$$f_3(4)=\min\{c_{4,6}+f_2(6);c_{2,7}+f_2(7)\}=18, P_3(4)=7;$$

digraph ex03 {
                rankdir=LR;
                size="8,5"
                node [shape = box];
                "1" -> "2" [ label = "3" ];
                "1" -> "3" [ label = "7" ];
                "1" -> "4" [ label = "2" ];
# Упростили     на 3 этапе
                "2" -> "6" [ label = "11" ];
                "3" -> "7" [ label = "7" ];
                "4" -> "7" [ label = "13" ];            
# Упростили     на 2 этапе
                "6" -> "8" [ label = "3" ];
                "7" -> "9" [ label = "1" ];
                "8" -> "10" [ label = "1" ];
                "9" -> "10" [ label = "4" ];
}
 

n=4

$$f_4(1)=\min\{c_{1,2}+f_3(2);c_{1,3}+f_3(3);c_{1,4}+f_3(4)\}=18, P_4(1)=3;$$

Таким образом оптимальный путь 1-2-6-8-10, затраты по которому составляют \(f_4(1)=18\)

digraph ex01 {
                rankdir=LR;
                size="8,5"
                node [shape = box];
                "1" -> "2" [ label = "3",style=bold,color=red ];
                "1" -> "3" [ label = "7",style=dotted];
                "1" -> "4" [ label = "2",style=dotted ];
                "2" -> "5" [ label = "9",style=dotted ];
                "2" -> "6" [ label = "11",style=bold,color=red ];
                "3" -> "5" [ label = "5",style=dotted ];
                "3" -> "6" [ label = "10",style=dotted ];
                "3" -> "7" [ label = "7",style=dotted ];
                "4" -> "6" [ label = "15",style=dotted ];
                "4" -> "7" [ label = "13",style=dotted ];
                "5" -> "8" [ label = "7",style=dotted ];
                "5" -> "9" [ label = "5",style=dotted ];
                "6" -> "8" [ label = "3",style=bold,color=red ];
                "6" -> "9" [ label = "4",style=dotted ];
                "7" -> "8" [ label = "7",style=dotted ];
                "7" -> "9" [ label = "1",style=dotted ];
                "8" -> "10" [ label = "1",style=bold,color=red ];
                "9" -> "10" [ label = "4",style=dotted ];
}

Решение в R

Для решения задач связанных с поиском кратчайшего пути в графе в GNU R можно использовать функцию пакет igraph. В таком случае

library(igraph)
mytable<-data.frame(from=c(1,1,1,2,2,3,3,3,4,4,5,5,6,6,7,7,8,9)
            to=c(2,3,4,5,6,5,6,7,6,7,8,9,8,9,8,9,10,10),
            weight=c(3,7,2,9,11,5,10,7,15,13,7,5,3,4,7,1,1,4))
## from to weight
## 1 1 2 3
## 2 1 3 7
## 3 1 4 2
## 4 2 5 9
## 5 2 6 11
## 6 3 5 5

Задаем граф, рисуем его и находим кратчайшие пути.

g<-graph.data.frame(mytable,directed =T )
plot(g)

distances(g,algorithm ="bellman-ford")

##     1  2  3  4  5  6  7  8  9 10
## 1   0  3  7  2 12 14 14 17 15 18
## 2   3  0 10  5  9 11 15 14 14 15
## 3   7 10  0  9  5 10  7 12  8 12
## 4   2  5  9  0 14 15 13 18 14 18
## 5  12  9  5 14  0  9  6  7  5  8
## 6  14 11 10 15  9  0  5  3  4  4
## 7  14 15  7 13  6  5  0  6  1  5
## 8  17 14 12 18  7  3  6  0  5  1
## 9  15 14  8 14  5  4  1  5  0  4
## 10 18 15 12 18  8  4  5  1  4  0

shortest_paths(g, 1, 10)$vpath

## [[1]]
## + 5/10 vertices, named:
## [1] 1  2  6  8  10

\[Z = 2500 х_1 + 3500 х_2 \to \max \]

\[\left\{ {\begin{array}{}
{3{x_1} + 10{x_2} \le 330}\\
{16{x_1} + 4{x_2} \le 400}\\
{6{x_1} + 6{x_2} \le 240}\\
{{x_1} \ge 0}\\
{{x_2} \ge 12}
\end{array}} \right.\]

Решение задачи начинается с построения области допустимых решений. При этом возможны следующие случаи:

1. Область допустимых решений — пустое множество. В этом случае решения нет из-за несовместимости ограничений.
2. Область допустимых решений  — единственная точка. Тогда она и является оптимальным решением.
3. Область допустимых решений — выпуклая неограниченная область. В таком случае решение может не существовать, если нет ограничений сверху при задаче на максимум, или снизу при задаче на минимум, а может находится в одной из угловых точек.
4. Область допустимых значений — выпуклый многоугольник. В этом случае можно найти координаты всех угловых точек, вычислить в них значение и выбрать оптимальное.

Однако существует иной способ. Пусть \(c_0\) — некоторое число. Прямая \(c_1x_1+c_2x_2=c_0\) является линией уровня целевой функции. В каждой точке этой прямой целевая функция принимает одно и то же значение, равное \(c_0\). Вектор — градиент целевой функции:

\[ \bar{c}= grad\;L(\bar{x})=\left(\frac{\partial L}{\partial x_1}; \frac{\partial L}{\partial x_2}\right)=(c_1,c_2) \]

перпендикулярен линиям уровня и показывает направление, в котором функция возрастает с наибольшей скоростью. Выбирая из линий уровня, проходящих через область допустимых значений, наиболее удаленную в направлении вектора \(\bar{c}\) (в случае минимизации — в противоположном направлении), определяем угловую точку, в котором целевая функция принимает максимальное (минимальное) значение. Если экстремум достигается сразу в двух смежных угловых точках, то по теореме об альтернативном оптимуме, оптимальным решением будет любая точка отрезка соединяющие эти точки.

Алгоритм графического метода

1. Построить область допустимых значений
2. Построить вектор градиент \(\bar{c}=(c_1,c_2) \)
3. Построить семейство линий уровня перпендикулярных вектору \(\bar{c}\), проходящую через область допустимых решений.
4. Выбрать линию уровня, проходящую через область допустимых решений наиболее удаленную в направлении вектора \(\bar{c}=(c_1,c_2) \)  (в задаче на максимум, при решении задачи на минимум — в противоположном направлении). Определить угловые точки области, через которые они проходят.
5. Найти координаты точек экстремума и значения целевой функции в этих точках

Решение задачи представлено на рисунке:

# Первый график

x1<- (-100:500)/10
old<-par(mar=c(1,1,1,1))
plot(0,type="n",xlab="",ylab="", xlim=c(-10, 50),ylim = c(-10, 50),bty="n",xaxt="n",yaxt="n")
polygon(c(0,0,10,20,22),c(12,33,30,20,12), col = "lightblue", border = NA)
axis(1,pos=c(0,0),at=c(-10,10,20,30,40))
axis(2,pos=c(0,0),las=2,at=c(-10,10,20,30,40))
arrows(-10.5,0,51,0,angle=15)
arrows(0,-10.5,0,51,angle=15)
lines(x1,(330-3*x1)/10,col="blue")
text(30,30,expression(3*х[1]+10*х[2]==330),cex=0.8,col="blue")
lines(x1,(400-16*x1)/4,col="red")
text(20,45,expression(16*х[1]+4*х[2]==400),cex=0.8,col="red")
lines(x1,(240-6*x1)/6,col="green")
text(8,40,expression(6*х[1]+6*х[2]==240),cex=0.8,col="green")
abline(h=12)
lines(x1,(-25*x1)/35,lty=3,lwd=3)
text(-5,5,expression(2500*х[1]+3500*х[2]==0),cex=0.8)
arrows(0,0,12.5,17.5)
lines(x1,(130000-2500*x1)/3500,lty=2,lwd=3)
text(35,20,expression(2500*х[1]+3500*х[2]==130000),cex=0.8)
points(10,30,cex=1.5,col="red",pch=19)

# Второй график

x1<- (-50:50)/10
old<-par(mar=c(1,1,1,1))
plot(0,type="n",xlab="",ylab="", xlim=c(-5, 5),ylim = c(-5, 5),bty="n",xaxt="n",yaxt="n")
polygon(c(-8,-8,8),c(9,-8,-7), col = "lightblue", border = NA)
axis(1,pos=c(0,0),at=c(-4,-3,-2,-1,1,2,3,4))
axis(2,pos=c(0,0),las=2,at=c(-4,-3,-2,-1,1,2,3,4))
arrows(-5.5,0,5.5,angle=15)
arrows(0,-5.5,0,5.5,angle=15)
abline(1,-1)
text(-3,2,expression(х[1]+х[2]<=1),col="black")
text(5,0.5,expression(х[1]),col="black")
text(0.5,5,expression(х[2]),col="black")

# Третий график

x1<- (-50:50)/10
old<-par(mar=c(1,1,1,1))
plot(0,type="n",xlab="",ylab="", xlim=c(-5, 5),ylim = c(-5, 5),bty="n",xaxt="n",yaxt="n")
polygon(c(-8,8,8),c(-7,9,-8), col = "lightblue", border = NA)
axis(1,pos=c(0,0),at=c(-4,-3,-2,-1,1,2,3,4))
axis(2,pos=c(0,0),las=2,at=c(-4,-3,-2,-1,1,2,3,4))
arrows(-5.5,0,5.5,angle=15)
arrows(0,-5.5,0,5.5,angle=15)
abline(1,1)
text(2,-3,expression(х[1]+х[2]<=1),col="black")
text(5,0.5,expression(х[1]),col="black")
text(0.5,5,expression(х[2]),col="black")

# Четвертый график
x1<- (-50:50)/10
old<-par(mar=c(1,1,1,1))
plot(0,type="n",xlab="",ylab="", xlim=c(-5, 5),ylim = c(-5, 5),bty="n",xaxt="n",yaxt="n")
polygon(c(-8,-8,8),c(-9,9,7), col = "lightblue", border = NA)
axis(1,pos=c(0,0),at=c(-4,-3,-2,-1,1,2,3,4))
axis(2,pos=c(0,0),las=2,at=c(-4,-3,-2,-1,1,2,3,4))
arrows(-5.5,0,5.5,angle=15)
arrows(0,-5.5,0,5.5,angle=15)
abline(-1,1)
text(-2,1,expression(х[1]-х[2]<=1),col="black")
text(5,0.5,expression(х[1]),col="black")
text(0.5,5,expression(х[2]),col="black")

Стандартная форма

В стандартной форме задаются \( n \) действительных чисел \(с_1, c_2, \ldots, c_n \); \( m \) действительных чисел \( b_1,b_2,\ldots ,b_m \); и \( mn \) действительных чисел \( a_{ij} \),  где \( i=1,2, \ldots, m \) и \(  j = 1,2, \ldots, n \).
Требуется найти \( n \) действительных чисел \( x_1,x_2,\ldots,x_n \) которые:

Максимизируют целевую функцию \( \sum\limits_{j=1}^n c_j x_j \) при заданных ограничениях:   \( \sum\limits_{j=1}^n a_{ij} x_j \le b_i\)  при \( i=1,2, \ldots, m \) и ограничениях неотрицательности \( x_j \ge 0 \) при \( j= 1,2,\ldots,n \)

Преобразование в стандартную форму

Задача находится не в стандартной форме если:

1. Целевая функция минимизируется, а не максимизируется
2. На имеющиеся переменные не наложены условия неотрицательности
3. Некоторые ограничения имеют форму равенства, т.е. имеют знак равенства вместо меньше или равно
4. Некоторые ограничения вместо знака «меньше или равно» имеют знак «больше или равно».

Как получить стандартную форму в таких случаях? Рассмотрим по пунктам:

1. Достаточно поменять знаки целевой функции, например: \( -5x_1+2x_2 \to min \) эквивалентно \( 5x_1-2x_2 \to max \).
2. Если отсутствует ограничение неотрицательности для переменной \(x_j\), то заменяем эту переменную выражением \( x_j^{‘}-x_j^{«} \) из двух неотрицательных переменных \(x_j^{‘},x_j^{«}  \ge 0 \).
3. Если ограничение имеет вид равенства, то заменяем его парой из двух неравеств «меньше или равно» и «больше или равно».
4. Для смены смены выда неравенства с «больше или равно» на «меньше или равно» умножаем обе части неравества на -1 и меняем знак сравнения.

Пример:

\[-5x_1+3x_2 \to min \]

\[ \left\{ {\begin{array}{} {x_1 + 2x_2 = 9 } \\ {x_1 — 3x_2 \le 7 } \\ {x_1 \ge 0 } \end{array}} \right. \]

Шаг 1. Меняем знак целевой функции

\[5x_1-3x_2 \to max \]

\[ \left\{ {\begin{array}{} {x_1 + 2x_2 = 9 } \\ {x_1 — 3x_2 \le 7 } \\ {x_1 \ge 0 } \end{array}} \right. \]

Шаг 2. Для второй переменной \(x_2\) ограничений неотрицательности нет. Заменим на выражение \(x_2=x_2^{‘}-x_2^{»} \):

\[5x_1-3x_2^{‘}+3x_2^{»} \to max \]

\[ \left\{ {\begin{array}{} {x_1 + 2x_2^{‘}-2x_2^{»} = 9 } \\ {x_1 — 3x_2^{‘}+3x_2^{»} \le 7 } \\ {x_1,x_2^{‘},x_2^{»} \ge 0 } \end{array}} \right. \]

Шаг 3. Заменяем равенство на два неравенства:

\[5x_1-3x_2^{‘}+3x_2^{»} \to max \]

\[ \left\{ {\begin{array}{} {x_1 + 2x_2^{‘}-2x_2^{»} \ge 9 } \\ {x_1 + 2x_2^{‘}-2x_2^{»} \le 9 }\\ {x_1 — 3x_2^{‘}+3x_2^{»} \le 7 } \\ {x_1,x_2^{‘},x_2^{»} \ge 0 } \end{array}} \right. \]

Шаг 4. Изменяем знак неравества:

\[5x_1-3x_2^{‘}+3x_2^{»} \to max \]

\[ \left\{ {\begin{array}{} {-x_1-2x_2^{‘}+2x_2^{»} \le -9 } \\ {x_1 + 2x_2^{‘}-2x_2^{»} \le 9 }\\ {x_1 — 3x_2^{‘}+3x_2^{»} \le 7 } \\ {x_1,x_2^{‘},x_2^{»} \ge 0 } \end{array}} \right. \]

Получили задачу в стандартном виде.

Примечание:

В данной заметке рассмотрен  основной подход к виду стандартной формы задачи ЛП. Существует и другой подход, например в книге Соколов А. В., Токарев В.В. Методы оптимальных решений. В 2-х т. Т 1. Общие положения. математическое программирование. — 3-е изд. — 564 с. на с. 423 выделяется стандартный вид задачи на максимум:

\[ \left\{ {\begin{array}{} {с_1 x_1 +\ldots+ c_n x_n \to max  } \\{a_{11} x_1 +\ldots+ a_{1n} x_n \le b_1 } \\ {\ldots} \\{a_{m1} x_1 +\ldots+ a_{mn} x_n \le b_m }\\{ x_1 ,\ldots ,x_n \ge 0 }\end{array}} \right. \]

и на минимум

\[ \left\{ {\begin{array}{} {с_1 x_1 +\ldots+ c_n x_n \to min  } \\{a_{11} x_1 +\ldots+ a_{1n} x_n \ge b_1 } \\ {\ldots} \\{a_{m1} x_1 +\ldots+ a_{mn} x_n \ge b_m }\\{ x_1 ,\ldots ,x_n \ge 0 }\end{array}} \right. \]

Однако в общепринят первый вариант из которого второй получается простым умножением обеих частей неравенств на -1.

Николай Кузнецов управляет небольшим механическим заводом. В будущем месяце он планирует изготавливать два продукта (А и В), по которым удельная маржинальная прибыль оценивается в 2500 и 3500 руб., соответственно.

Изготовление обоих продуктов требует затрат на машинную обработку, сырье и труд На изготовление каждой единицы продукта А отводится 3 часа машинной обработки, 16 единиц сырья и 6 единиц труда. Соответствующие требования к единице продукта В составляют 10, 4 и 6. Николай прогнозирует, что в следующем месяце он может предоставить 330 часов машинной обработки, 400 единиц сырья и 240 единиц труда. Технология производственного процесса такова, что не менее 12 единиц продукта В необходимо изготавливать в каждый конкретный месяц.

Николай хочет построить модель с тем, чтобы определить количество единиц продуктов А и В, которые он доложен производить в следующем месяце для максимизации маржинальной прибыли.

Решение задачи

Этап 1. Определение переменных

Существует целевая переменная (обозначим её Z), которую необходимо оптимизировать, то есть максимизировать или минимизировать (например, прибыль, выручка или расходы). Николай стремится максимизировать маржинальную прибыль, следовательно, целевая переменная:

Z — это суммарная маржинальная прибыль (в рублях), полученная в следующем месяце в результате производства продуктов А и В.

Существует ряд неизвестных искомых переменных (обозначим их х₁, х₂, х₃ и пр.), чьи значения необходимо определить для получения оптимальной величины целевой функции, которая, в нашем случае является суммарной маржинальной прибылью. Эта маржинальная прибыль зависит от количества произведенных продуктов А и В. Значения этих величин необходимо рассчитать, и поэтому они представляют собой искомые переменные в модели. Итак, обозначим:

х₁ — количество единиц продукта А, произведенных в следующем месяце.

х₂ — количество единиц продукта В, произведенных в следующем месяце.

Очень важно четко определить все переменные величины; особое внимание уделите единицам измерения и периоду времени, к которому относятся переменные.

Этап. 2. Построение целевой функции

library(lpSolve) # Подключили библиотеку
f.obj <- c(2500, 3500) # Описали целевую функцию
names(f.obj) <-c("A","B")
a.mat<-rbind(c(3,10), # матрица
             c(16,4),   # коээфициентов
             c(6,6),    # при ограничениях
             c(1,0),
             c(0,1))   
a.dir<-c("<=","<=","<=",">=",">=")
b.vec<-c(330,400,240,0,12) # вектор ограничений

result<-lp ("max", f.obj, a.mat, a.dir, b.vec)

result

## Success: the objective function is 130000

result$solution

## [1] 10 30

library(linprog)

## Warning: package 'linprog' was built under R version 3.2.2

(result<-solveLP( f.obj, b.vec, a.mat, TRUE,const.dir=a.dir,lpSolve=T))

## 
## 
## Results of Linear Programming / Linear Optimization
## (using lpSolve)
## 
## Objective function (Maximum): 130000 
## 
## Solution
##   opt
## A  10
## B  30
## 
## Constraints
##   actual dir bvec free
## 1    330  <=  330    0
## 2    280  <=  400  120
## 3    240  <=  240    0
## 4     10  >=    0   10
## 5     30  >=   12   18

(mtx<-matrix(c(0,1,2,-1,0,1,-2,1,0),ncol=3))

##      [,1] [,2] [,3]
## [1,]    0   -1   -2
## [2,]    1    0    1
## [3,]    2    1    0

min(mtx[1,])

## [1] -2

min(mtx[2,])

## [1] 0

min(mtx[3,])

## [1] 0

max(min(mtx[1,]),min(mtx[2,]),min(mtx[3,]))

## [1] 0

max(mtx[,1])

## [1] 2

max(mtx[,2])

## [1] 1

max(mtx[,3])

## [1] 1

min(max(mtx[,1]),max(mtx[,2]),max(mtx[,3]))

## [1] 1

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    8    6    4    4    3
## [2,]    5    3    2    2    1
## [3,]    4    7    7    3    5
## [4,]    5    3    2    2    1
## [5,]    1    4    4    2    3

(a<-apply(mtx,1,min))

## [1] 3 1 3 1 1

max(a)

## [1] 3

(b<-apply(mtx,2,max))

## [1] 8 7 7 4 5

min(b)

## [1] 4

mtx

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    8    6    4    4    3
## [2,]    5    3    2    2    1
## [3,]    4    7    7    3    5
## [4,]    5    3    2    2    1
## [5,]    1    4    4    2    3

mtx[c(1,3),]

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    8    6    4    4    3
## [2,]    4    7    7    3    5

(mtx1<-mtx[c(1,3),c(3,4,5)])

##      [,1] [,2] [,3]
## [1,]    4    4    3
## [2,]    7    3    5

max(apply(mtx1,1,min))

## [1] 3

min(apply(mtx1,2,max))

## [1] 4

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    3    7    1    1    5
## [2,]    4    9    3    6    2
## [3,]    2    3    1    4    7

(a<-max(apply(mtx,1,min)))

## [1] 2

(b<-min(apply(mtx,2,max)))

## [1] 3

library(lpSolve)
(result<-lp("min",c(1,1,1), t(mtx), rep(">=",5),c(1,1,1,1,1)))

## Success: the objective function is 0.3684211

result$objval

## [1] 0.3684211

result$solution

## [1] 0.00000000 0.31578947 0.05263158

(a<-result$solution/result$objval)

## [1] 0.0000000 0.8571429 0.1428571

(result<-lp("max",c(1,1,1,1,1), mtx, rep("<=",3),c(1,1,1)))

## Success: the objective function is 0.3684211

result$objval

## [1] 0.3684211

result$solution

## [1] 0.0000000 0.0000000 0.2631579 0.0000000 0.1052632

(b<-result$solution/result$objval)

## [1] 0.0000000 0.0000000 0.7142857 0.0000000 0.2857143

# Функция, возвращающая индекс стратегии
get.k<-function(vec){
  cusum=0
  tst=runif(1)
  for(i in 1:length(vec)) {
   if(vec[i]==0) next
   cusum<-cusum+vec[i]
   if(tst>cusum) next
  return(i)
  }
}
set.seed(2015)
test<-c()
for(i in 1:1000) test=c(test,mtx[get.k(a),get.k(b)])

table(test)

## test
##   1   2   3   7 
##  86 221 656  37

summary(test)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   2.000   3.000   2.755   3.000   7.000

test<-c()
for(i in 1:1000) test=c(test,mtx[sample(1:3, 1),get.k(b)])

table(test)

## test
##   1   2   3   5   7 
## 478  79 249  97  97

summary(test)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   1.000   2.000   2.547   3.000   7.000

test<-c()
for(i in 1:1000) test=c(test,mtx[get.k(a),sample(1:5, 1)])

table(test)

## test
##   1   2   3   4   6   7   9 
##  22 183 211 203 164  22 195

summary(test)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   3.000   4.000   4.726   6.000   9.000

Решение в LibreOffice

Пример решения приведен в файле matrix.ods.

Решение сводится к решениям прямой и обратной задач линейного программирования с помощью встроенной системы Решатель.

Результаты получаются аналогичными вышеприведенному решению.

Рассмотрим задачу. Рекламное агентство, в штате которого десять человек получило заказ на рекламу нового продукта на радио и телевидении. Данные о рекламной аудитории, стоимости рекламы и количестве занятых при её изготовлении агентов заданы в таблице.

Сколько минут рекламного времени должно купить агентство на радио и ТВ, чтобы максимизировать аудитоию и минимизировать издержки, если контракт запрещает более 6 минут на радио?

Имеем следующую задачу:

$$\left\{ {\begin{array} {} {{u_1} = 4{x_1} + 8{x_2} \to \max }\\ {{u_2} = 8{x_1} + 24{x_2} \to \min }\\ {{x_1} \le 6}\\ {{x_1} + 2{x_2} \le 10}\\ {{x_1} \ge 0;\;{x_2} \ge 0} \end{array}} \right.$$

Если решать задачу только на максимум, то имеем классическую задачу линейного программирования, которая решается следующим образом (предполагается, что библиотека lpSolve установлена).

library(lpSolve) 
f.obj <- c(4, 8) # Описали целевую функцию
names(f.obj) <-c("X1","X2")
a.mat<-rbind(c(1,0), # матрица
             c(1,2),   # коээфициентов
             c(1,0),    # при ограничениях
             c(0,1))   
a.dir<-c("<=","<=",">=",">=")
b.vec<-c(6,10,0,0) # вектор ограничений
(result<-lp ("max", f.obj, a.mat, a.dir, b.vec))

## Success: the objective function is 40

result$solution

## [1] 0 5

Если решать эту же задачу только на минимум:

library(lpSolve) 
f.obj <- c(4, 24) # Описали целевую функцию
names(f.obj) <- c("X1","X2")
a.mat <- rbind(c(1,0), # матрица
             c(1,2),   # коээфициентов
             c(1,0),    # при ограничениях
             c(0,1))   
a.dir <- c("<=","<=",">=",">=")
b.vec <- c(6,10,0,0) # вектор ограничений
(result <- lp ("min", f.obj, a.mat, a.dir, b.vec))

## Success: the objective function is 0

result$solution

## [1] 0 0

Графически задача может быть представлена следующим образом

x1 <- (-10:100)/10
old <- par(mar=c(1,1,1,1))
plot(0,type="n",xlab="",ylab="", xlim=c(-1, 10),ylim = c(-1, 10),bty="n",xaxt="n",yaxt="n")
grid()
polygon(c(0,0,6,6),c(0,5,2,0), col = "lightblue", border = NA)
axis(1,pos=c(0,0),at=c(-1,1,2,3,4,5,6,7,8,9))
axis(2,pos=c(0,0),las=2,at=c(-1,1,2,3,4,5,6,7,8,9))
arrows(-1.2,0,10.1,0,angle=15)
arrows(0,-1.2,0,10.1,angle=15)
lines(x1,(10-x1)/2,col="blue")
text(4,4,expression(х[1]+2*х[2]==10),cex=0.8,col="blue")
abline(v=6)
text(0.5,10,expression(х[2]))
text(10,-0.5,expression(х[1]))
points(0,0,cex=1.5,col="red",pch=19)
points(0,5,cex=1.5,col="red",pch=19)

par(old)

Множество всех решений может быть представлено:

x1<-seq(0,6,by=0.1)
x2<-seq(0,6,by=0.1)
d<-expand.grid(x1=x1,x2=x2)
d<-subset(d,x1+2*x2<=10)
d$u1<-4*d$x1+8*d$x2
d$u2<-8*d$x1+24*d$x2
plot(d$u2~d$u1,type="p",pch=19,col="blue",main="Множество решений",xlab=expression(u[1]),ylab=expression(u[2]),cex=0.8)
points(0,0,cex=1.5,col="red",pch=19)
points(40,120,cex=1.5,col="red",pch=19)
points(20,50,cex=1,col="yellow",pch=19)
text(21,55,"A",col="yellow")
points(23,53,cex=1,col="yellow",pch=19)
text(24,55,"B",col="yellow")

points(20,40,cex=1,col="green",pch=19)
text(21,36,"C")

Где u₁ — аудитория в миллионах человек (эффективность), а u₂ — стоимость рекламы. Рассмотрим решения А (20,50) и B(23,53). Вариант A имеет меньшую стоимость , но вариант B более эффективен. В таком случае можно говорить, что варианты несравнимы. Рассмотрим точку C(20,40). Как видим, при той же аудитории стоимость данного решения меньше.

Найдем множество всех таких точек и построим на графике.

plot(d$u2~d$u1,type="p",pch=19,col="lightblue",main="Множество решений",xlab=expression(u[1]),ylab=expression(u[2]),cex=0.8)
z <- aggregate(u1~u2,data=d,max)
z <- aggregate(u2~u1,data=z,min)
lines(z$u2~z$u1,col="blue",lwd=2)

Данная линия называется Парето-оптимальными вариантами.

(u1<-max(d$u1))

## [1] 40

(u2<-min(d$u2))

## [1] 0

plot(d$u2~d$u1,type="p",pch=19,col="lightblue",main="Множество решений",xlab=expression(u[1]),ylab=expression(u[2]),cex=0.8)
points(u1,u2,pch=19,col="red")

d.new<-d
l<-2
d.new$s<-((u1-d.new$u1)^l+(u2-d.new$u2)^l)^(1/l)
(u.opt<-d.new[which.min( d.new$s),])

##    x1 x2 u1 u2        s
## 21  2  0  8 16 35.77709

plot(d$u2~d$u1,type="p",pch=19,col="lightblue",main="Множество решений",xlab=expression(u[1]),ylab=expression(u[2]),cex=0.8)
points(u1,u2,pch=19,col="red")
arrows(u1,u2,u.opt$u1,u.opt$u2)
text(u.opt$u1+2,u.opt$u2+5,paste0("U(",u.opt$u1,";",u.opt$u2,")"))

(z<-max(d$u1))

## [1] 40

head(d.new<-subset(d,u1==z))

##       x1  x2 u1    u2
## 1281 6.0 2.0 40  96.0
## 1340 5.8 2.1 40  96.8
## 1399 5.6 2.2 40  97.6
## 1458 5.4 2.3 40  98.4
## 1517 5.2 2.4 40  99.2
## 1576 5.0 2.5 40 100.0

(u.opt<-d.new[which.min( d.new$u2),])

##      x1 x2 u1 u2
## 1281  6  2 40 96

plot(d$u2~d$u1,type="p",pch=19,col="lightblue",main="Множество решений",xlab=expression(u[1]),ylab=expression(u[2]),cex=0.8)
lines(d.new$u1,d.new$u2,col="blue",lwd=2)
points(u.opt$u1,u.opt$u2,pch=19,col="red",cex=1.5)
text(u.opt$u1-3,u.opt$u2-5,paste0("U(",u.opt$u1,";",u.opt$u2,")"))

d.new<-d
d.new$s<-d.new$u1*0.8-d.new$u2*0.2
head(d.new)

##    x1 x2  u1  u2    s
## 1 0.0  0 0.0 0.0 0.00
## 2 0.1  0 0.4 0.8 0.16
## 3 0.2  0 0.8 1.6 0.32
## 4 0.3  0 1.2 2.4 0.48
## 5 0.4  0 1.6 3.2 0.64
## 6 0.5  0 2.0 4.0 0.80

tail(d.new)

##       x1  x2   u1    u2    s
## 2931 0.2 4.8 39.2 116.8 8.00
## 2932 0.3 4.8 39.6 117.6 8.16
## 2990 0.0 4.9 39.2 117.6 7.84
## 2991 0.1 4.9 39.6 118.4 8.00
## 2992 0.2 4.9 40.0 119.2 8.16
## 3051 0.0 5.0 40.0 120.0 8.00

subset(d.new,s==max(d.new$s))

##      x1 x2 u1 u2    s
## 1281  6  2 40 96 12.8

d.new<-d
d.new$s<-d.new$u1*0.5-d.new$u2*0.5
head(d.new)

##    x1 x2  u1  u2    s
## 1 0.0  0 0.0 0.0  0.0
## 2 0.1  0 0.4 0.8 -0.2
## 3 0.2  0 0.8 1.6 -0.4
## 4 0.3  0 1.2 2.4 -0.6
## 5 0.4  0 1.6 3.2 -0.8
## 6 0.5  0 2.0 4.0 -1.0

tail(d.new)

##       x1  x2   u1    u2     s
## 2931 0.2 4.8 39.2 116.8 -38.8
## 2932 0.3 4.8 39.6 117.6 -39.0
## 2990 0.0 4.9 39.2 117.6 -39.2
## 2991 0.1 4.9 39.6 118.4 -39.4
## 2992 0.2 4.9 40.0 119.2 -39.6
## 3051 0.0 5.0 40.0 120.0 -40.0

subset(d.new,s==max(d.new$s))

##   x1 x2 u1 u2 s
## 1  0  0  0  0 0

k1<-0.8
k2<-0.2
f.obj <- c(4*k1-8*k2, 8*k1-24*k2) # Описали целевую функцию
names(f.obj) <-c("X1","X2")
a.mat<-rbind(c(1,0), # матрица
             c(1,2),   # коээфициентов
             c(1,0),    # при ограничениях
             c(0,1))   
a.dir<-c("<=","<=",">=",">=")
b.vec<-c(6,10,0,0) # вектор ограничений
(result<-lp ("max", f.obj, a.mat, a.dir, b.vec))

## Success: the objective function is 12.8

result$solution

## [1] 6 2

Для работы с приведенными примерами необходимо, чтобы был установлен R, а также RStudio. Ссылки на ПО есть на соответствующей странице.

Для работы необходимо установить некоторые пакеты (просто скопируйте приведенные команды в консоль и нажмите Enter)

install.packages("ahp")
install.packages("shinythemes")
install.packages("shinyAce")
install.packages("shinyjs")

После установки пакетов вы готовы к запуску графической оболочки. Наберите команды.

library(ahp)
RunGUI()

Если все сделано правильно, должно запуститься окно браузера.

Может возникнуть ошибка, вызванная тем, что какой-то из пакетов не встал, например, из-за сетевого сбоя. Если обратить внимание на текст сообщения при установке, то можно заметить, что пакет ahp при инсталляции также ставит необходимые для своей работы пакеты: jsonlite, httpuv, xtable, htmlwidgets, shiny, rstudioapi, visNetwork, data.tree, formattable, DiagrammeR. Если вместо окна браузера появляется ошибка, то надо просто установить недостающий пакет командой install.packages  или из меню приложения.

Рассмотрим решение задачи выбора лидера из статьи, описанной в Wikipedia.

В данной задаче необходимо выбрать из трех кандидатов одного на должность руководителя. Кандидаты оцениваются по критериям: возраст, опыт, образование и личные качества. На рисунке показана иерархия для этой задачи. Простейшая иерархия содержит три уровня: цель, критерии и альтернативы. Числа на рисунке показывают приоритеты элементов иерархии с точки зрения цели, которые вычисляются в МАИ на основе парных сравнений элементов каждого уровня относительно связанных с ними элементами вышерасположенного уровня. Приоритеты альтернатив относительно цели (глобальные приоритеты) вычисляются на заключительном этапе метода путем линейной свертки локальных приоритетов всех элементов. В данном примере лучшим кандидатом является Дик, так как имеет максимальное значение глобального приоритета.

Данный пример уже имеется в системе, достаточно выбрать его из выпадающего меню (tom_dick_harry.ahp)

Перейдем на вкладку Visualize и посмотрим на схему задачи

На вкладке анализ мы увидим результаты расчетов

Аналогичные результаты можно получить и из командной строки.

ahpFile <- system.file("extdata", "tom_dick_harry.ahp", package="ahp")
tomAhp <- Load(ahpFile)
Visualize(tomAhp)

Calculate(tomAhp) # Пересчитать
Analyze(tomAhp)  # Результаты в текстовом виде
# Weight  Dick   Tom Harry Inconsistency
# 1 Choose the Most Suitable Leader 100.0% 48.1% 38.5% 13.4%          4.4%
# 2  ¦--Experience                   54.8% 34.9% 14.1%  5.7%          3.3%
# 3  ¦--Charisma                     27.0%  5.2% 20.1%  1.7%          6.1%
# 4  ¦--Education                    12.7%  4.2%  2.8%  5.6%            NA
# 5  °--Age                           5.6%  3.8%  1.5%  0.4%          2.5%
AnalyzeTable(tomAhp) # Результаты в графическом виде

Известные проблемы и пути их решения

Главная проблема — пока не удалось заставить работать с русским языком в Windows. В Mac OS X все работает прекрасно, a в Windows скрипт вылетает по ошибке. Проблема понятна и она заключается в кодовых страницах, поскольку R понимает latin1 и UTF-8, а в OS X кодировка одна, UTF-8, то проблем нет. В Windows же кодировка  windows-1251 и иногда возникают ошибки, в случае, если приложение не адаптировано к иным языкам, отличным от английского.

Поэтому, как вариант, можно использовать обходное решение. Называть все элементы латиницей. а потом в документ уже вставлять рисунки на русском языке. Для построения дерева программа использует graphviz и можно получить описание дерева с помощью вызова функции.

cat(GetGraph(tomAhp)$dot_code)

Получаем описание схемы на языке dot, которую можно отредактировать простыми поиском и заменой в любом текстовом редакторе

digraph {

  'Choose the Most Suitable Leader' [label = 'Choose the Most Suitable Leader', style = 'filled,rounded', shape = 'box', color = 'honeydew4', fillcolor = 'honeydew', penwidth = '4', fontname = 'helvetica'] 
  'Experience' [label = 'Experience', style = 'filled,rounded', shape = 'box', color = 'honeydew4', fillcolor = 'honeydew', penwidth = '4', fontname = 'helvetica'] 
  'Education' [label = 'Education', style = 'filled,rounded', shape = 'box', color = 'honeydew4', fillcolor = 'honeydew', penwidth = '4', fontname = 'helvetica'] 
  'Charisma' [label = 'Charisma', style = 'filled,rounded', shape = 'box', color = 'honeydew4', fillcolor = 'honeydew', penwidth = '4', fontname = 'helvetica'] 
  'Age' [label = 'Age', style = 'filled,rounded', shape = 'box', color = 'honeydew4', fillcolor = 'honeydew', penwidth = '4', fontname = 'helvetica'] 
  'Alternatives' [label = 'Alternatives', style = 'filled,rounded', shape = 'box', color = 'thistle4', fillcolor = 'thistle', penwidth = '4', fontname = 'helvetica'] 
  'Tom' [tooltip = 'age: 50
education: 4
experience: 7
leadership: 10
', style = 'filled,rounded', shape = 'box', color = 'thistle4', fillcolor = 'thistle', penwidth = '4', fontname = 'helvetica'] 
  'Dick' [tooltip = 'age: 60
education: 6
experience: 10
leadership: 6
', style = 'filled,rounded', shape = 'box', color = 'thistle4', fillcolor = 'thistle', penwidth = '4', fontname = 'helvetica'] 
  'Harry' [tooltip = 'age: 30
education: 8
experience: 5
leadership: 6
', style = 'filled,rounded', shape = 'box', color = 'thistle4', fillcolor = 'thistle', penwidth = '4', fontname = 'helvetica'] 
'Choose the Most Suitable Leader'->'Experience' [arrowhead = 'vee', color = 'grey35', penwidth = '2'] 
'Choose the Most Suitable Leader'->'Education' [arrowhead = 'vee', color = 'grey35', penwidth = '2'] 
'Choose the Most Suitable Leader'->'Charisma' [arrowhead = 'vee', color = 'grey35', penwidth = '2'] 
'Choose the Most Suitable Leader'->'Age' [arrowhead = 'vee', color = 'grey35', penwidth = '2'] 
'Experience'->'Alternatives' [arrowhead = 'vee', color = 'grey35', penwidth = '2'] 
'Education'->'Alternatives' [arrowhead = 'vee', color = 'grey35', penwidth = '2'] 
'Charisma'->'Alternatives' [arrowhead = 'vee', color = 'grey35', penwidth = '2'] 
'Age'->'Alternatives' [arrowhead = 'vee', color = 'grey35', penwidth = '2'] 
'Alternatives'->'Tom' [color = 'grey35', penwidth = '2', dir = 'none'] 
'Alternatives'->'Dick' [color = 'grey35', penwidth = '2', dir = 'none'] 
'Alternatives'->'Harry' [color = 'grey35', penwidth = '2', dir = 'none'] 
}

@Manual{ahp2016,
  title = {ahp: Analytic Hierarchy Process},
  author = {Christoph Glur},
  year = {2016},
  note = {R package version 0.2.8},
  url = {https://CRAN.R-project.org/package=ahp}
}
@Book{Saati1993,
  Title                    = {Принятие решений. Метод анализа иерархий},
  Address                  = {М.},
  Author                   = {Саати, Т.},
  Language                 = {russian},
  Pages                    = {278},
  Numpages                 = {278},
  Publisher                = {Радио и связь},
  Year                     = {1993}
}