rplacroix 16 июл в 10:32

Алгоритм Тарьяна для поиска минимального набора уравнений

Средний

6 мин

1.7K

Дана система, состоящая из большого количества уравнений (необязательно линейных), где вам необходимо найти всего лишь несколько переменных. Как это сделать эффективно? Какой минимальный набор уравнений вам потребуется? В этой статье мы обсудим графовое представление систем уравнений, применим алгоритм Тарьяна и формализуем процесс на Python.

Постановка задачи

Введём двудольный граф $G = (V_1 \cup V_2, \ E)$ . V_1 — множество вершин, представляющие уравнения, V_2 — вершины, представляющие переменные. Ребро $e \in E$ соединяет пар, если уравнение зависит от переменной .

И будем пользоваться матрицей смежности . По строкам — уравнения, по стобцам - переменные. Если $M_{ij} = 1$ , то -ое уравнение зависит от -переменной.

$S:\left\{\begin{array}{cccc} f_1\left(x_1, x_2, \ldots, x_n\right) & = & 0 & e_1 \\ f_2\left(x_1, x_2, \ldots, x_n\right) & = & 0 & e_2 \\ \vdots & & & \\ f_n\left(x_1, x_2, \ldots, x_n\right) & = & 0 & e_n \end{array},\right.$

Найти минимальную подсистему $v_1 \subset V_1$ уравнений, решив которую, можно найти переменные $v_2 = \{x_i, x_j, \dots\} \subset V_2$ (произвольное подмножество).

Формулировка алгоритма

Количество уравнений и переменных должно совпадать для применения данного алгоритма.

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

$M^{'} = \begin{pmatrix} 1 & a_{12} & \dots & a_{1n}\\ a_{21} & 1 & \dots & a_{2n} \\ \vdots & \vdots & \dots & a_{n - 1, n} \\ \dots & \dots & \dots & 1 \end{pmatrix}$

На самом деле, задача состоит в приведении матрицы к блочной нижнетреугольной (БНТ) форме. Это мы сделаем с помощью алгоритма Тарьяна, который найдёт компоненты сильной связности (КСС) в ориентированном графе. КСС - это такой подграф, что для любой пары вершин А и Б, существует путь из А в Б и обратно - из Б в А.

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

Hidden text

def tarjan_algorithm_with_start(graph: dict[int, set[int]], start_node: int):
    def dfs(v):
        nonlocal index
        indices[v] = index
        lowlink[v] = index
        index += 1
        stack.append(v)
        on_stack[v] = True

        for neighbor in graph[v]:
            if indices[neighbor] == -1:
                dfs(neighbor)
                lowlink[v] = min(lowlink[v], lowlink[neighbor])
            elif on_stack[neighbor]:
                lowlink[v] = min(lowlink[v], indices[neighbor])

        if indices[v] == lowlink[v]:
            scc = []
            while True:
                node = stack.pop()
                on_stack[node] = False
                scc.append(node)
                if node == v:
                    break
            strongly_connected_components.append(scc)

    index = 0
    stack = []
    indices = [-1] * len(graph)
    lowlink = [-1] * len(graph)
    on_stack = [False] * len(graph)
    strongly_connected_components = []

    dfs(start_node)

    return strongly_connected_components

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

C точки зрения системы уравнений, добавляется вспомогательная переменная $r = f_r(x_1, \dots x_h)$ , выраженная некой функцией от интересующих нас переменных.

Пример

$\left\{\begin{array}{l}\dot{x_1}(t)=-x_1(t) \\\dot{x_2}(t)=x_1(t)-x_2(t) \\\dot{x_3}(t)=x_1(t) \\y_1(t)=3 x_2(t)+x_1(t) \\y_2(t)=2 x_3(t)\end{array} \right.$

Три переменных состояния ( x_1, x_2, x_3 ) и две алгебраические переменные y_1, y_2 . Введём вспомогательных уравнения, показывающие связь между переменной x_i и её производной $\dot{x_i}$ :

$\left\{\begin{array}{lll}\dot{x_1}+x_1 & =0 & e_0 \\\dot{x_2}-x_1+x_2 & =0 & e_1 \\\dot{x_3}-x_1 & =0 & e_2 \\y_1-3 x_2-x_1 & =0 & e_3 \\y_2-2 x_3 & =0 & e_4 \\x_1-f\left(\dot{x_1}\right) & =0 & e_5 \\x_2-g\left(\dot{x_2}\right) & =0 & e_6 \\x_3-h\left(\dot{x_3}\right) & =0 & e_7\end{array}\right.$

Матрица смежности:

$\begin{aligned}\begin{array}{l} e_0 \\ e_1 \\ e_2 \\ e_3 \\ e_4 \\ e_5 \\ e_6 \\ e_7 \end{array}\left(\begin{array}{llllllll} 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \end{array}\right) \\ & \end{aligned}$

Переставляем строки так, чтобы на главной диагонали были 1. Программа (см. ниже раздел Код, где представлена наивная реализация такого алгоритма) даст перестановку $\{0, 1, 7, 5, 6, 2, 3, 4\}$ :

$\begin{aligned}& \begin{array}{l}e_0 \\e_1 \\e_7 \\e_5 \\e_6 \\e_2 \\e_3 \\e_4\end{array}\left(\begin{array}{llllllll}1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\1 & 1 & 0 & 0 & 0 & 0 & 1 & 0 \\0 & 0 & 1 & 0 & 0 & 0 & 0 & 1\end{array}\right) \\&\end{aligned}$

В системе восемь переменных, перечислим их здесь для удобства:

$\begin{array}{llllllllll} x_1 & x_2 & x_3 & \dot{x}_1 & \dot{x_2} & \dot{x}_3 & y_1 & y_2 \end{array}$

Предположим, что нужно найти минимальный набор уравнений, которые нужны для вывода переменной y_1 . Вывод программы: [[3, 0], [4, 1], [6]].

Рис. 1. Граф с выделенными КСС, необходимых для нахождения

Уравнения дадут нужный результат. Здесь мы подействовали на номера [[3, 0], [4, 1], [6]]перестановкой {0, 1, 7, 5, 6, 2, 3, 4}: 3 переходит в 5, 0 переходит в 0, 4 переходит в 6 и т.д.

Записав уравнения в заданном порядке получим матрицу сниженнной размерности в блочной нижнетреугольной форме: (переменные / столбцы в таком порядке: $x_1, \dot{x_1}, x_2, \dot{x_2}, y_1$ ).

$\begin{aligned}& \begin{array}{l}e_5 \\e_0 \\e_6 \\e_1 \\e_3\end{array}\end{aligned}\left(\begin{array}{lllll}1 & 1 & 0 & 0 & 0 \\1 & 1 & 0 & 0 & 0 \\0 & 0 & 1 & 1 & 0 \\0 & 1 & 1 & 1 & 0 \\0 & 1 & 0 & 1 & 1\end{array}\right)$

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

$\left\{\begin{array}{l}\dot{x_1}(t)=-x_1(t) \\\dot{x_2}(t)=x_1(t)-x_2(t) \\y_1(t)=3 x_2(t)+x_1(t)\end{array} \right.$

Пример 2

Найдём уравнения, необходимые для вычисления y_1 и y_2 . Введём корневую вершину (на рисунке обозначена 8):

Рис. 2. Граф с добавленной корневой вершиной

Вывод программы: [[3, 0], [4, 1], [6], [5], [2], [7], [8]]. Потребуются все уравнения в системе, чтобы найти . Однако алгоритм привёл матрицу в удобную БНТ форму.

При подготовке данного материала использовался пример и концепция из статьи Minimal Equation Sets for Output Computation in Object-Oriented Models. Vincenzo Manzoni Francesco Casella Dipartimento di Elettronica e Informazione, Politecnico di Milano.

Код

Hidden text

import itertools
import numpy as np
from numpy.typing import NDArray

from tarjan.utils import tarjan_algorithm_with_start, matr2dct, colorGraph


class MinimalSubsetDetector:
    def __init__(self, matrix: NDArray, omitOnes: bool = True):
        assert matrix.shape[0] == matrix.shape[1], 'Матрица должна быть квадратной'
        permuted = MinimalSubsetDetector.permuteMatrix(matrix)
        if permuted is None:
            raise ValueError('У переданной матрицы не существует такой перестановки строк, что на главной диагонали стоят единицы.')
        self.matrix = permuted
        self.size = self.matrix.shape[0]
        self.hasRoot = False
        if omitOnes:
            for i in range(self.size):
                self.matrix[i][i] = 0

    @staticmethod
    def permuteMatrix(matrix) -> NDArray | None:
        """
        Находит такую перестановку строк, что на главной диагонали должны стоять только 1.
        O(k^N), где k - примерное кол-во 1 в строке.
        :param matrix:
        :return:
        """
        result = []
        for row in matrix:
            indices = np.where(row == 1)[0]
            result.append(indices.tolist())
        result_dict = {}
        for i, sublist in enumerate(result):
            for element in sublist:
                if element not in result_dict:
                    result_dict[element] = [i]
                else:
                    result_dict[element].append(i)
        result_list = [result_dict.get(i, []) for i in range(len(result))]
        combinations = itertools.product(*result_list)
        valid_perm = next(comb for comb in combinations if len(set(comb)) == len(comb))
        if valid_perm is None:
            print('Такой перестановки не существует.')
            return
        print(f'Необходимая перестановка: {valid_perm}.')
        return matrix[list(valid_perm)]

    @classmethod
    def initfromTeX(cls, dimension: int, texcode: str, omitOnes: bool = True):
        """
        :param omitOnes: опустить единицы на главной диагонали матрицы (убрать петли).
        :param dimension: размерность матрицы
        :param texcode: код для матрицы
        :return:
        """
        tex_matrix = texcode.replace("\\", '&').replace("\n", "")
        elements = tex_matrix.split("&")
        int_elements = [int(element.strip()) for element in elements if len(element)]
        return cls(np.array(int_elements).reshape(dimension, dimension), omitOnes=omitOnes)

    def addRootNode(self, interested_nodes: list[int]):
        """
        Добавление корневой вершины R, которая связана ребрами с interested_nodes. Применяется, когда
        нужно узнать минимальный набор уравнений для нахождения группы неизвестных.
        :param interested_nodes:
        :return:
        """
        self.hasRoot = True
        row = [0] * self.size
        for n in interested_nodes:
            row[n] = 1
        new_row = np.array(row)
        extended_matrix = np.vstack([self.matrix, new_row])
        new_column = np.array([0] * (self.size + 1))
        new_column = new_column[:, np.newaxis]
        self.matrix = np.hstack((extended_matrix, new_column))

    def find(self, node: int = None):
        if node is None and not self.hasRoot:
            raise ValueError('Укажите индекс переменной, которую вы хотите найти.')
        return tarjan_algorithm_with_start(matr2dct(self.matrix), start_node=self.size if self.hasRoot else node)

    def color_SCC(self, scc):
        return colorGraph(self.matrix, scc)


def test1():
    detector = MinimalSubsetDetector.initfromTeX(
        8,
        r'1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\'
        r'1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\'
        r'0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\'
        r'1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\'
        r'0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\'
        r'1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\'
        r'1 & 1 & 0 & 0 & 0 & 0 & 1 & 0 \\'
        r'0 & 0 & 1 & 0 & 0 & 0 & 0 & 1'
    )
    # Для поиска минимального набора уравнений для переменных у1 и у2.
    # detector.addRootNode([6, 7])
    answers = detector.find(6)  # найти для у1
    print(answers)
    detector.color_SCC(answers)


if __name__ == '__main__':
    test1()

Утилитные функции:

Hidden text

import networkx as nx
import numpy as np
from matplotlib import pyplot as plt


def matr2dct(m):
    gr = {}
    num_nodes = m.shape[0]
    for i in range(num_nodes):
        neighbors = set(np.nonzero(m[i])[0])
        gr[i] = neighbors
    return gr


def colorGraph(matrix, scc):
    G = nx.DiGraph(matrix)
    pos = nx.spring_layout(G)
    colors = ['red', 'green', 'blue', 'yellow', 'orange', 'violet', 'pink']
    for component, color in zip(scc, colors):
        nx.draw_networkx_nodes(G, pos, nodelist=component, node_color=color, node_size=500)
    nx.draw_networkx_edges(G, pos, arrows=True)
    nx.draw_networkx_labels(G, pos)
    plt.axis('off')
    # plt.savefig('graph.png')
    plt.show()

Теги:

Хабы:

Алгоритм Тарьяна для поиска минимального набора уравнений

Постановка задачи

Формулировка алгоритма

Пример

Пример 2

Код

Публикации

Истории

Работа

Ближайшие события