9. Графы и сети

9.0. Введение

Сети - математические объекты, обладающие дискретной и нерегулярной структурой. Структура задается графом сети - совокупностью узлов и ветвей. Ветви могут быть однонаправленными (ориентированный граф) или двунаправленными (неориентированный граф).

Сеть отличается от графа тем, что с ее ветвями (а иногда и узлами) связаны некоторые численные характеристики. Сети используются при моделировании разветвленных коммуникаций, механизмов, технологических процессов, процессов принятия решений.

В этой главе приводится несколько примеров, демонстрирующих, как можно программировать на "МК-85" небольшие сетевые задачи.

9.1. Связывающая сеть

Программа 144. Алгоритм Прима получает из полносвязной сети с двунаправленными ветвями связную сеть минимальной стоимости.

Входные данные: n - число узлов сети; c_ij, (j>i) - стоимость ветви (непосредственной связи) между узлами i и j.

Выходные данные: {i_k, j_k}, k=1...n-1 - индексы ветвей, обеспечивающих минимум суммарной стоимости сети при сохранении ее связности.

Ограничение размерности: n≤16.

10 INPUT A:DEFM (A-1)*A/2-17
20 FOR B=1 TO A-1:FOR C=B+1 TO A:PRINT B;C;
30 INPUT I((C-2)*(C-1)/2+B):NEXT C:NEXT B:D=1
40 I=1E20:FOR B=1 TO A:IF FRAC (D/2-B)<.5 THEN 90
50 FOR C=1 TO A:F=0:IF C<B;F=(B-2)*(B-1)/2+C
60 IF C>B;F=(C-2)*(C-1)/2+B
70 IF I(F)<I;I=I(F):G=B:H=C:E=F
80 NEXT C
90 NEXT B:I(E)=1E20
100 IF FRAC (D/2-H)<.5;D=D+2-H/2:PRINT G;H;I
110 IF D<2-A-1 THEN 40

Размер: 265, I((n-1)*n/2)

Пример: n=6; C=((10, 24, 66, 91, 15), (..., 17, 22, 51, 43), (..., 52, 37, 29), (...,21, 66), (..., 81)).

9.2. Геодезическая сеть

Программа 145. Уравновешивание невязок углов в свободной геодезической сети. (Другие задачи уравновешивания требуют модификации программы лишь в части вычисления невязок.)

Коррелаты невязок полигонов xi определяются путем решения системы уравнений:

где k - число полигонов;
n_i - число точек в i-том полигоне;
n_ij - число смежных точек полигонов i и j;
u_li - отсчет внутреннего угла (в градусах) на l-той точке i-го полигона.

Все отсчеты сети - равноточные. Невязки углов можно будет выразить через коррелаты полигонов по следующим формулам:

для внешней точки i-го полигона - Du=x_i;
для смежной точки полигонов i и j - Du=x_i-x_j;
для узловой точки - Du=x_i-0,5еx_j, j¦i.

(Эти формулы не включены в программу. )

При вводе данных предполагается, что точки сети пронумерованы и их общее число равно m. Ввод данных осуществляется в порядке обхода точек. Для каждой точки указываются номера смежных полигонов и соответствующие внутренние углы. Если точка внутренняя, последний угол может вычисляться автоматически дополнением до 360. Для перехода к новой точке вводите i=0 (если точка внутренняя, можете вводить u=0).

Изменения точности вычислений можно добиться путем изменения порога чувствительности в строке 150.

10 VAC:INPUT "m",A,"k",B:FOR C=1 TO A:D=360:PRINT C;
20 FOR E=1 TO B:G(E)=0:NEXT E
30 INPUT "i",E:IF E=0 THEN 60
40 G(E)=1:E=E+B:INPUT "u",F:IF F>0;D=D-F:G(E)=G(E)+F:GOTO 30
50 G(E)=G(E)+D
60 FOR E=1 TO B:FOR F=E TO B:IF G(E)+G(F)<2 THEN 80
70 GOSUB 180:G(G)=G(G)+1
80 NEXT F:NEXT E:NEXT C:FOR E=1 TO B:FOR F=E TO B:GOSUB 180
90 IF E¦F;G(G)=1-G(G):GOTO 110
100 G(E+B)=(G(G)-2)*180-G(E+B):G(E)=0
110 NEXT F:NEXT E
120 $="":FOR C=1 TO B:A=G(C+B):FOR D=1 TO B:IF C=D THEN 150
130 E=C:F=D:IF F<E;E=D:F=C
140 GOSUB 180:A=A-G(G)*G(D)
150 NEXT D:F=C:GOSUB 180:A=A/G(G):IF ABS (G(C)-A)>1E^-10;$=CHR 1
160 G(C)=A:NEXT C:IF $¦"" THEN 120
170 FOR A=1 TO B:PRINT A;G(A):NEXT A:END
180 G=(F-1)*F/2+E+B+B:RETURN

Размер: 475, G((k+5)*k/2)

Пример:

Ввод данных:

Ответ: x=(-4,11764705; -0,352941175; -0,6470588258)^T

9.3. Сеть систем массового обслуживания

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

Рассмотрим конкретный пример:

Число систем в сети n=3. Узел 0 (S)- внешний источник заявок; узлы 1, 2, 3 - системы массового обслуживания. p_ij вероятности перехода заявки из системы i в систему j.

Программа 146. Оценка средних времен ожидания.

Интенсивности потоков заявок на входах систем

Коэффициенты передачи узлов a_i определяются из системы уравнений:

Допустимая размерность 2≤n≤9.

Узел сети - многоканальная система массового обслуживания с неограниченным общим буфером. Загрузка такой системы

где J_i - среднее время обслуживания;
m_i - число каналов обслуживания.

Если r_i оказывается больше единицы, программа запрашивает другие J_i и m_i.

Среднее время ожидания:

DEFM n*(n+1)-19
10 INPUT "n",A:FOR B=0 TO A:FOR C=1 TO A:PRINT B;C;
20 INPUT D:IF B=0;D=-D
30 IF C=B;D=D-1
40 G(A*B+C)=D:NEXT C:NEXT B
50 FOR B=1 TO A-1:FOR C=B+1 TO A
60 F=G(B*A+B):G=G(B*A+C):IF F=0;IF G=0 THEN 100
70 E=SQR (F*F+G*G):D=G/E:E=F/E
80 FOR F=0 TO A:G=E*G(F*A+B)+D*G(F*A+C)
90 G(F*A+C)=E*G(F*A+C)-D*G(F*A+B):G(F*A+B)=G:NEXT F
100 NEXT C:NEXT B:FOR B=A TO 1 STEP -1:D=0:IF B=A THEN 120
110 FOR C=B+1 TO A:D=D+G(A*A+C)*G(C*A+B):NEXT C
120 G(A*A+B)=(G(B)-D)/G(B*A+B):NEXT B
130 INPUT "l",C:FOR B=1 TO A:H=G(A*A+B)*C:PRINT B;"l";H
140 INPUT "y",E,"m",F:E=E*H:PRINT "r";E/F:IF E≥F THEN 140
150 I=0:G=1/E:D=0
160 D=D+G:I=I+1:IF I<F;G=G*(F-I)/E:GOTO 160
170 D=E/H/(F-E+(F-E)-2*D):PRINT "w";D:NEXT B:GOTO 130

Размер: 508, G(n*n+n)

Расчет примера.

Ввод: n=3; l₀=0,5; P= (1; 0; 0); (0,1; 0,3; 0,6); (0; 0; 0,2); (0; 0; 0))
Вывод: l₁=0,5556.
Ввод: J₁=10; m₁=6.
Вывод: r₁=0,9259; w₁=18,09; l₂=0,1667.
Ввод: J₂=7; m₂=2.
Вывод: r₂=0,5833; w₂=3,611; l₃=0,3667.
Ввод: J₃=1; m₃=1.
Вывод: r₃=0,3667; w₃=0,5789.

9.4. Граф вероятных состояний

Моделируемая динамическая система может находиться в любом из случайных состояний S_i, i=1...n. Известны интенсивности l_ij переходов из состояния i в состояние j, а также исходное состояние S_k (p_k(0)=1). Необходимо проследить вероятности состояний p_i(t) во времени.

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

Программа 147. Система уравнений Колмогорова. В данном случае для решения используется метод Рунге - Кутта 4-го порядка (см. описание программы 115). Помимо данных, перечисленных выше, необходимо задать шаг интегрирования s. Для завершения списка переходов вводите i=0. В процессе счета программа высвечивает модельное время, не прерываясь. Чтобы вывести на дисплей текущий вектор вероятностей, достаточно нажать любую клавишу (за исключением [EXE]) и немного подождать.

DEFM 4*n+3*m-20
10 VAC:INPUT "n",A:A=A*4:B=A
20 INPUT "i",D:IF D≤0 THEN 40
30 B=B+3:D(B)=D*4:INPUT "j",D,"l",F(B):E(B)=D*4:GOTO 20
40 INPUT "k",C,"s",F:C(4*C)=1
50 E=E+F:PRINT E;:FOR C=4 TO A STEP 4:E(C)=C(C):F(C)=0:NEXT C
60 FOR D=0 TO 1.5 STEP .5:FOR C=4 TO A STEP 4:D(C)=0:NEXT C
70 FOR C=A+3 TO B STEP 3:D(D(C))=D(D(C))-F(C)*C(D(C)):$=$+KEY
80 D(E(C))=D(E(C))+F(C)*C(D(C)):NEXT C:FOR C=4 TO A STEP 4
90 F(C)=F(C)+F*D(C)*(1+SGN FRAC (D/1.5))
100 C(C)=C(C)+F*D(C)*(1+INT D)/2:NEXT C:NEXT D:CSR 0
110FOR C=4 TO A STEP 4:C(C)=E(C)+F(C)/6:IF $¦"";PRINT C/4;C(C)
120 NEXT C:$="":GOTO 50

Размер: 498, F(4*A+3*B)

Пример:

Ввод данных:

DEFM 4*4+3*6-20
SET 7
4
1 2 3.4
2 3 7.2
3 4 2.6
4 1 1.3
1 3 2.1
3 1 3.4
0
1 0.1

Ответ: В момент t=0,2 p=(0,5296883; 0,1313437; 0,2429235; 0,09604452)^T

9.5. Герт-сеть

Ветви Герт-сети обозначают действия, характеризуемые матожиданием m_i и дисперсией D_i времени выполнения. Узлы отображают ветвление процесса, описываемое распределением вероятностей перехода p_i. На выходе каждого узла еp_i=1.

Программа 148. Правило Мейсона для спектральных коэффициентов передачи ветвей:

где W_t - коэффициент передачи всей сети;
E - множество контуров, содержащих фиктивную обратную связь;
L - множество контуров, не содержащих фиктивную обратную связь; (фиктивная обратная связь - это условно вводимая ветвь, соединяющая сток сети с ее истоком);
i - номер ветви;
W_i - коэффициент передачи ветви;
C - контур (кольцо ветвей);
k(C) - порядок контура; (контур порядка k - это k непересекающихся простых контуров).

Подстановка в уравнение Мейсона вероятностей p_i позволяет вычислить вероятность p_t перехода от источника к стоку, а двукратное аналитическое дифференцирование уравнения позволяет записать аналогичные выражения для матожидания m_t=f{p_i,m_i} и дисперсии D_t=f{p_i,m_i,D_i}. Эти три выражения заложены в программу.

Программа хранит в памяти список ветвей:

Ограничение размерности: n≤41.

Выделять контуры будем визуально (в противном случае программа оказалась бы слишком длинной). Список контуров зависит от того, между какими узлами определяется коэффициент передачи. Между выбранной парой узлов следует замкнуть фиктивную обратную связь. При вводе списка: в ответ на запрос программы "k" указывается порядок контура; в ответ на запрос "e" контуры, содержащие фиктивную ветвь, помечаются цифрой "1", не содержащие - цифрой "0"; затем перечисляются номера ветвей, составляющих контур. Для перехода к следующему контуру вводите i=0, для окончательного завершения ввода - k=0.

10 INPUT "n",A:DEFM A*3-14:FOR B=1 TO A:C=B*3:PRINT B;
20 INPUT "p",J(C),"m",K(C),"D",L(C):NEXT B
30 C=1:D=0:E=0:F=0:G=0:H=0
40 J=1:K=0:L=0:INPUT "k",I:IF I=0 THEN 100
50 INPUT "e",B
60 INPUT "i",A:A=A*3:IF A=0 THEN 80
70 J=J*J(A):K=K+K(A):L=L+L(A):GOTO 60
80 I=-1-I:L=(L+K*K)*J:K=K*J:FOR A=0 TO 2
90 C(A+A+B)=C(A+A+B)+I*J(A):NEXT A:GOTO 40
100 J=-D/C:K=D*E/C/C-F/C:L=D*E/C-F:L=G*D-2*E*L:L=L/C/C-H/C
110 K=K/J:PRINT "p";J,"m";K,"D";L/J-K*K:GOTO 30

Размер: 354, L(3*C)

Пример:

Ввод списка ветвей:

n=5
p₁=1; m₁=3; D₁=6
p₂=0,2; m₂=1; D₂=1
p₃=0,1; m₃=2; D₃=3
p₄=0,4; m₄=0; D₄=0
p₅=0,3; m₅=4; D₅=0

1) Эквивалентная ветвь s - t:

2) Эквивалентная ветвь s - q:

9.6. Перт-сеть

Ветви Перт-сети - это действия, которые должны выполняться в заданной последовательности, включая и параллельное выполнение. Узел сети считается выполненным, если выполнены все входящие в него ветви. Задержки в ветвях описываются матожиданием m_i и дисперсией s_i². Алгоритм определения времени выполнения сети по методу Перт осуществляет суммирование каждой из оценок по путям от истока к стоку. В узлах, где пути сходятся, выбирается максимальная из сумм. Следует иметь в виду, что операция вычисления максимума m и s² не имеет строгого математического обоснования, и при наличии в сети параллельных ветвей с близкими характеристиками метод заметно занижает результирующую оценку матожидания и завышает дисперсию.

Программа 149. Оценка раннего времени выполнения. Узлы должны быть пронумерованы подряд, начиная с единицы, так чтобы все ветви были направлены от младших номеров к старшим. Количество узлов не должно превышать 83. Для задания нового узла ответьте "+" на запрос "E" (программа сообщит текущий номер узла n); для завершения ввода сети ответьте E=".". Для текущего узла перечисляются все заходящие в него ветви. На запрос "r" введите номер узла, из которого исходит эта ветвь, далее введите m_i и s_i. Для перехода к новому узлу введите r=0. Программа сообщит ожидемый момент выполнения текущего узла m_n и его СКО s_n.

DEFM 2k-21 (k - количество узлов в сети)
10 VAC:B=2
20 B=B+2:PRINT "n";B/2;
30 INPUT " r",A:IF A=0 THEN 70
40 INPUT "m",C,"s",D:C=C+C(A+A):IF C>C(B);C(B)=C
50 D=D*D+D(A+A):IF D>D(B);D(B)=D
60 GOTO 30
70 PRINT "m";C(B),"s";SQR D(B):INPUT "E",$:IF $¦"." THEN 20

Размер: 160, D(2*B)

Пример:

Список ветвей

Ввод-вывод данных:

9.7. Дерево решений

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

Программа 150. Задача об инвестициях. Адаптированная версия.

Вы - владелец коммерческой фирмы. Через надежные каналы до Вас дошла информация о том, что через неделю все излишки со счетов предприятий изымут в централизованный фонд помощи перестройке. У Вас горят R рублей. Вы озадачили друзей - и начался мозговой штурм.

Когда шторм утих, Вы изучили его дары.

В общей сложности друзья, а также приятели их знакомых предложили Вам n выгодных способов пожарного вложения денег, каждый из которых требует ровно s_i рублей и обещает ровно через год c_i дохода в твердой валюте, i=1...n. Но получился небольшой перебор - оказалось, что ваших "деревянных" не хватает на реализацию всех предложений. Приходится отбирать то, что по зубам. Разумеется, любой законопослушный гражданин будет сортировать не друзей, а предложения. Таким образом, задача может быть строго формализована:

Найти такое I={i} - подмножество {1...n}, для которого

Решение такой задачи "в лоб" требует перебора всех комбинаций. Предположим, что n=30. Тогда всего комбинаций 2³⁰, то есть миллиард - и у Вас мало шансов получить ответ при жизни. Методы дискретной оптимизации дают рецепт того, как уменьшить число просматриваемых вариантов, не упустив оптимальный. Так, не следует просматривать нереализуемые и заведомо неперспективные направления поиска.

Предварительно расположим предложения в порядке уменьшения их рентабельности - так чтобы выполнялось:

Такое упорядочивание обеспечивает скорейшее нахождение оптимума. Скорость зависит еще и от того, насколько точна оценка перспективности направления поиска. Здесь мы примем самую близорукую оценку - текущую сумму потерь дохода. В начале поиска примем потери субоптимального варианта большими (P=10⁹⁹), а текущие потери нулевыми (p₀=0).

Процесс поиска представляет собой чередование прямого и обратного хода. Прямой ход - i - i+1, обратный - i - i-1. Начинаем с i=1 и движемся вперед.

На прямом ходе будем принимать/отклонять предложения последовательно, начиная с первого. Мы можем принять текущее предложение, если на него хватает денег, то есть r_i-1≥s_i. В этом случае остаток средств уменьшится:

В противном случае будем отказываться от предложения. Но тогда увеличатся накопленные потери:

Если накопленные потери сравнялись с тем, что мы имеем для уже найденного претендента на лучший вариант,

дальнейшее движение вперед в этом направлении бесперспективно - нужно возвращаться. Но этого может и не произойти. Тогда, рассмотрев последнее, n-е предложение, мы записываем новый рекорд P и запоминаем тот выбор, который к нему привел. Оптимальный выбор записывается в виде бинарного вектора x. x_i=1, если i-е предложение было принято, и x_i=0, если отклонено.

Разумеется, двигаясь вперед, мы запоминали текущий выбор в бинарном векторе t. Таким образом, запоминание нового рекорда сводится к следующим присвоениям:

Как вести себя, двигаясь обратно?

Если текущий выбор был t_i=0, продолжаем двигаться обратно.

Если же t_i=1, возможно, что p_i<P. Если это действительно так, переключаемся на t_i=0 и пытаемся двигаться вперед. Если же, двигаясь назад, мы вернулись в i=0, поиск закончен.

Ограничение памяти - n≤22. Распределение переменных:

Для пересчета задачи с новым R можете использовать команду RUN 30.

10 INPUT A:DEFM 6*A-22:FOR B=1 TO A:D=B*6-6:PRINT "c";B;
20 INPUT F(D),"s",E(D):NEXT B:A=A*6-6
30 INPUT "R",C:B=1E99:D=-6
40 D=D+6:IF D>A THEN 100
50 IF E(D)≤C(D);G(D)=1:I(D)=C(D)-E(D):J(D)=D(D):GOTO 40
60 G(D)=0:I(D)=C(D):J(D)=D(D)+F(D):IF J(D)<B THEN 40
70 D=D-6:IF D<0 THEN 110
80 IF G(D)=1 THEN 60
90 GOTO 70
100 B=D(D):FOR D=0 TO A STEP 6:H(D)=G(D):NEXT D:GOTO 70
110 B=0:FOR D=0 TO A STEP 6
120 IF H(D)>0;B=B+F(D):C=C-E(D):PRINT D/6+1;
130 NEXT D:PRINT CHR 0,"C";B,"r";C

Размер: 339, D(n*6)

Пример: n=10;