Известия РАН. Теория и системы управления, 2019, № 4, стр. 62-73

СИНТЕЗ БЫСТРЫХ И СВЕРХБЫСТРЫХ РЕШАТЕЛЕЙ БОЛЬШИХ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ МЕТОДАМИ ТЕОРИИ УПРАВЛЕНИЯ

М. Г. Гаджиев¹, К. В. Жгун¹, Н. Е. Зубов^2, *, В. Н. Рябченко^1, 2

¹ МЭИ (национальный исследовательский ун-т)
Москва, Россия

² МГТУ им. Н.Э. Баумана
Москва, Россия

^* E-mail: nik.zubov@gmail.com

Поступила в редакцию 06.04.2018
После доработки 06.11.2018
Принята к публикации 26.11.2018

DOI: 10.1134/S0002338819020082

Полный текст (PDF)

Аннотация

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

0. Введение. Одним из основных объектов линейной алгебры являются большие системы линейных алгебраических уравнений (СЛАУ)

$Ax = b,$

где $A \in {{F}^{{n \times n}}}$ и $b \in {{F}^{n}}$ – заданные матрица и вектор соответственно. Здесь и в дальнейшем символ F обозначает либо поле вещественных чисел $\mathbb{R}$, либо поле комплексных чисел $\mathbb{C}$.

Известно большое число итерационных методов и программных пакетов итерационных решателей (iterative solvers), предназначенных для решения (0.1). Тем не менее следует отметить, что эти методы работают при достаточно жестких ограничениях на свойства матрицы А, среди которых необходимо отметить цикличность, нормальность, симметричность, положительную определенность [1–5]. Более того, до сих пор не известен итерационный решатель уравнений, который успешно справлялся с рассматриваемой задачей для общего случая матрицы А. Так методы, основанные на подпространствах Крылова, например, метод сопряженных градиентов или итераций порядка r GMRES(r) (A Generalized Minimal Residual Method), демонстрирует хорошую глобальную сходимость для положительно определенных и нормальных матриц [1].

В случае отсутствия свойства нормальности у матрицы А динамика итерационных процессов решения уравнения (0.1) и соответствующие условия сходимости методов на основе подпространств Крылова становятся достаточно сложными, в том числе, хаотическими и до конца не изученными. Так, в [6] представлены примеры простых уравнений (0.1), при решении которых итерации GMRES(1) быстро сходятся, достигая точности 10^–15, а ошибка процесса GMRES(2) демонстрирует хаотическое поведение в пределах 10⁰–10^–1.

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

Цель настоящей статьи – построение на основе оригинального метода многошаговой декомпозиции и синтеза управления линейной многомерной динамической системой [7] алгоритмов быстрых и сверхбыстрых итерационных решателей СЛАУ (0.1), обеспечивающих в численных расчетах быстро сходящиеся итерационные процессы. Другими словами, требуется, чтобы данные итерационные решатели обладали максимальным запасом устойчивости, а их сходимость наступала бы в течение первых итераций независимо от свойств матриц, представляющих СЛАУ. Минимизация сложности алгоритмов (complexity) в данном случае не рассматривается.

1. Подход к анализу итерационных решателей на основе методов теории управления. Параметры итерационных решателей, такие, как сдвиг, релаксация и т.д., можно рассматривать как управляющие воздействия, а сами решатели – управляемыми динамическими системами. Один из первых подходов в этом направлении, как указано в [1], можно проследить в ранней работе Р. Беллмана [8]: “…Вычисления можно рассматривать как процесс управления, в котором мы хотим каким-то соответствующим образом сочетать сложность, время и точность. Кроме того, этот процесс следует рассматривать во многих случаях как адаптивное управление, в котором результаты предыдущего расчета используются для последующих расчетов, обеспечивая выбор не только параметров, но и алгоритмов”² ².

В последние полтора десятилетия за рубежом был изучен ряд численных алгоритмов с позиций стабилизации динамических систем на основе синтеза обратной связи. В [1] к этому направлению исследований отнесены работы [9, 10], где для улучшения решателей обыкновенных дифференциальных уравнений (ODE-solvers) синтезировалось управляющее воздействие, обеспечивающее выбор подходящего шага интегрирования. В публикациях [11–14] понятие управляемости использовалось в отношении методов вычисления собственных значений матричных пучков и решателей линейных матричных уравнений (0.1). В [15–17] для оптимизации решателя СЛАУ применялись методы теории робастной устойчивости.

Кратко рассмотрим некоторые наиболее часто используемые итерационные решатели линейного матричного уравнения (0.1), построенные на основе метода простой итерации (метода Ричардсона), методов полиномиальной и линейной итерации [1].

Для обратимой матрицы A уравнение (0.1) имеет хорошо известное решение

(1.1)

$x = {{A}^{{ - 1}}}b.$

Итерационный решатель на основе метода Ричардсона можно отнести к дискретной билинейной управляемой системе с одним входом и многими выходами (SIMO – Single Input Multi Outputs), заданной над Fⁿ:

(1.2)

${{x}_{{k + 1}}} = {{x}_{k}} + {{u}_{k}}(b - A{{x}_{k}}),\quad {{x}_{0}} \in {{F}^{n}}.$

Здесь ${{u}_{k}}$ – управляющий параметр (скалярное управление), ${{x}_{0}}$ – начальное условие. Неподвижная (установившаяся) точка в n-мерном пространстве этой итерации является решением (1.1) линейного уравнения (0.1).

Предложены различные стратегии для ${{u}_{k}}$ и некоторых семейств матриц [18–21]. В частности, при ${{u}_{k}} = u = {\text{const}}$ дискретная система (1.2) имеет установившееся значение (эквивалентно, является асимптотически устойчивой), если множество собственных значений, обозначаемое как ${\text{eig}}$, матрицы ${{E}_{n}} - uA$, где ${{E}_{n}}$ – единичная матрица порядка n, лежит внутри единичного круга на комплексной плоскости, т.е.

(1.3)

$\forall \lambda \in {\text{eig(}}{{E}_{n}} - uA{\text{)}},\quad {\text{|}}\lambda {\text{|}} < 1.$

Если же ${{u}_{k}}$ соответствует закону обратной связи, то

(1.4)

${{u}_{k}} = \frac{{e_{k}^{{\text{т }}}A{{e}_{k}}}}{{{\text{||}}A{{e}_{k}}{\text{|}}{{{\text{|}}}^{2}}}},$

где

(1.5)

${{e}_{k}} = b - A{{x}_{k}}$

– текущая невязка. Такой подход используется в методе на основе подпространств Крылова GMRES(1). При этом известно, что динамическая система (1.2) с управлением (1.4) асимптотически устойчива, а итерационный процесс сходится, если матрица $A + {{A}^{{\text{т }}}}$ положительно определена. Очевидно, что это условие не достаточно для произвольных матриц.

Обобщением метода Ричардсона на основе полиномиальных итераций порядка m является решатель в виде следующей дискретной системы:

(1.6)

${{x}_{{k + 1}}} = ({{E}_{n}} - {{p}_{k}}(A)A){{x}_{k}} + {{p}_{k}}(A)b,\quad {{x}_{0}} \in {{F}^{n}}.$

Здесь скалярное управление ${{p}_{k}}(A)$ есть полином степени не более m, например полином Чебышева. В последнем случае получается известный алгоритм m-циклического метода Ричардсона с чебышевскими параметрами [22]. Исследования свойств сходимости системы (1.6) и выбор ${{p}_{k}}(A)$ выполнены в работах [6, 23–25].

В [17] введена альтернатива билинейной SIMO-системе (1.6) в виде линейной системы с многими входами и многими выходами (MIMO – Multi Inputs Multi Outputs) следующего вида:

(1.7)

${{x}_{{k + 1}}} = ({{E}_{n}} - A){{x}_{k}} + G{{u}_{k}} + b,\quad {{x}_{0}} \in {{F}^{n}},$

где $G$ – некоторая заданная матрица, в общем случае не имеющая прямого отношения к уравнению (0.1), ${{u}_{k}}$ – вектор управления. Асимптотическую устойчивость системы и соответственно сходимость итерационного процесса (1.7) предлагалось обеспечивать с помощью закона обратной связи

(1.8)

${{u}_{k}} = K{{x}_{k}},$

синтезированного на основе линейно-квадратической оптимизации. Это в свою очередь приводит к известному алгоритму решателя LQRES [1].

2. Управляемость итерационных решателей. Известно, что необходимым и достаточным условием успешности синтеза закона управления с обратной связью вида (1.8) выступает условие полной управляемости динамической системы. Поэтому естественным в контексте настоящей работы является анализ управляемости рассматриваемых итерационных решателей.

В [1] выполнен анализ управляемости итерационного решателя, построенного на основе метода Ричардсона (1.2). Показано, что для любой входной последовательности $({{u}_{k}})$ траектория состояния $({{x}_{k}})$ сходится к решению уравнения (0.1), если и только если последовательность невязок

(2.1)

${{e}_{{k + 1}}} = ({{E}_{n}} - {{u}_{k}}A){{e}_{k}}$

сходится к нулю. Таким образом, динамика дискретной системы (1.2) эквивалентна динамики (2.1).

Справедливы следующие утверждения [1].

Утверждение 1. Существует открытое множество обратимых циклических вещественных (заданных над $\mathbb{R}$) матриц A размера $n \times n$, для которых итерационный решатель (1.2) на основе метод Ричардсона полностью управляем. В частности, это происходит, если A имеет n различных вещественных собственных значений.

Утверждение 2. Итерационный решатель (1.2) на основе метод Ричардсона не полностью управляем, если матрица A имеет одно чисто мнимое собственное значение.

Утверждение 3. Если вещественная матрица A имеет чисто мнимые собственные значения и ${{x}_{0}} \in {{\mathbb{R}}^{n}}$, то матричное уравнение (0.1) может быть решено любым неполиномиальным итерационным решателем на основе метода Ричардсона.

Напомним, что открытое множество – это множество, в котором каждый элемент входит в него вместе с некоторой окрестностью (в метрических пространствах и, в частности, на числовой прямой).

Далее, полиномиальный решатель (1.6) степени m ≥ 2 является управляемым для обратимых циклических матриц A [1].

Управляемость итерационного решателя в виде дискретной MIMO-системы (1.7) определяется с помощью известных критериев Калмана, Попова-Белевича-Хотиса и/или критериев на основе ленточных матриц [7]. Ясно, что (1.7) полностью управляема тогда и только тогда, когда полностью управляемой является пара матриц $({{I}_{n}} - A{\text{,}}G)$. Например, в терминах рангового критерия Калмана это соответствует равенству

(2.2)

${\text{rank[}}G({{E}_{n}} - A)G \cdots {{({{E}_{n}} - A)}^{{n - 1}}}G{\text{]}} = n.$

3. Линейно-квадратическое управление итерационными решателями. Другой подход к разработке итерационных методов решения линейного матричного уравнения (0.1) основан на методах синтеза линейных регуляторов. Ключевыми здесь являются следующие рассуждения [1].

Пусть заданы вещественные матрица $A \in {{\mathbb{R}}^{{n \times n}}}$ и вектор $b \in {{\mathbb{R}}^{n}}$ и требуется найти матрицу $G \in {{\mathbb{R}}^{{n \times m}}}$ и глобально убывающую последовательность векторов (управление) ${{u}_{k}} \in {{\mathbb{R}}^{m}}$, $k \in \mathbb{N}$ (здесь $\mathbb{N}$ – множество натуральных чисел), т.е.

$\mathop {\lim }\limits_{k \to \infty } {{u}_{k}} = 0,$

а также, что итерации на основе MIMO-системы (1.7) сходятся к решению (1.1).

Без потери общности предположим $b \in \operatorname{Im} G$ (вектор принадлежит образу преобразования с матрицей G, другими словами, уравнение $b = Gz$ разрешимо относительно вектора z). В противном случае матрица G может быть заменена на матрицу $\tilde {G} = [b\,\,G] \in {{\mathbb{R}}^{{n \times m + 1}}}$, для которой, как видно, выполняется последнее условие.

В предположении обратимости матрицы A решение (1.1) может быть записано в известном виде степенного ряда [5]

$x = {{A}^{{ - 1}}}b = \sum\limits_{j = 0}^{n - 1} {{{\alpha }_{j}}{{{({{E}_{n}} - A)}}^{j}}b} $

для некоторых ${{\alpha }_{j}} \in \mathbb{R}$, $j = 0,1, \ldots ,n - 1$. Таким образом, обратная матрица ${{A}^{{ - 1}}}$ может быть определена с помощью MIMO-системы (1.7). При этом поведение невязки (1.5) описывается дискретной линейной системой [1]

(3.1)

${{e}_{{k + 1}}} = ({{E}_{n}} - A){{e}_{k}} - AG{{u}_{k}},\quad {{e}_{0}} \in {{\mathbb{R}}^{n}}.$

На основании анализа (3.1) можно сделать следующие утверждения.

Утверждение 4. Дискретная MIMO-система (3.1) полностью управляемая тогда и только тогда, когда полностью управляемой является пара $({{E}_{n}} - A{\text{,}} - {\kern 1pt} AG)$ в смысле критерия (2.2).

Утверждение 5. В случае полной управляемости MIMO-системы (3.1) существует (при m > 1 неединственный) закон обратной связи (1.8), что выполняется условие асимптотической устойчивости:

(3.2)

$\forall \lambda \in {\text{eig(}}{{E}_{n}} - A + AGK{\text{)}},\quad {\text{|}}\lambda {\text{|}} < 1.$

Существует несколько подходов к формированию закона (1.8). Так, в [1] предлагается синтезировать подходящую в смысле (3.2) матрицу K на основе теории линейно-квадратичного гауссовского (ЛКГ) регулятора. В этом случае указанная матрица имеет вид

(3.3)

$K = - {{({{E}_{n}} + {{G}^{{\text{т }}}}{{A}^{{\text{т }}}}PAG)}^{{ - 1}}}{{G}^{{\text{т }}}}{{A}^{{\text{т }}}}P({{E}_{n}} - A).$

Здесь P – единственное симметрическое положительно определенное решение алгебраического матричного уравнения Риккати следующего вида:

(3.4)

$P = {{E}_{n}} + {{\tilde {A}}^{{\text{т }}}}P\tilde {A} + {{({{\tilde {G}}^{{\text{т }}}}P\tilde {A})}^{{\text{т }}}}{{({{E}_{m}} + {{\tilde {G}}^{{\text{т }}}}P\tilde {G})}^{{ - 1}}}{{\tilde {G}}^{{\text{т }}}}P\tilde {A},$

где $\tilde {A} = {{E}_{n}} - A$, $\tilde {G} = - AG$.

Тогда ЛКГ-алгоритм итерационного решателя LQRES [1] принимает вид:

Шаг 1. Выбрать такую матрицу G, что пара $({{E}_{n}} - A{\text{,}} - \,AG)$ полностью управляема.

Шаг 2. Вычислить решение P матричного уравнения Риккати (3.4).

Шаг 3. Вычислить матрицу K согласно выражению (3.3).

Шаг 4. Сформировать итерационный решатель вида

(3.5)

${{x}_{{k + 1}}} = ({{E}_{n}} - A){{x}_{k}} + GK(b - A{{x}_{k}}) + b.$

В соответствии с ЛКГ-теорией имеет место следующее соотношение для невязки и управления:

(3.6)

$\sum\limits_{k = 0}^\infty {({\text{||}}{{e}_{k}}{\text{|}}{{{\text{|}}}^{2}} + \,\,{\text{||}}{{u}_{k}}{\text{|}}{{{\text{|}}}^{2}})} = e_{0}^{{\text{т }}}P{{e}_{0}}.$

Отсюда вытекает, что последовательность невязок ${{e}_{k}}$ сходится к нулю, а траектория состояния ${{x}_{k}}$ – к решению (1.1).

В [1] показаны недостатки и преимущества LQRES. Во-первых, хотя выбор матрицы G на шаге 1 и не является тривиальным для произвольной матрицы A, он всегда может быть выполнен. Во-вторых, именно подходящая матрица G обеспечивает повышение скорости сходимости итерационного процесса. В-третьих, требование управляемости пары $({{E}_{n}} - A{\text{,}} - {\kern 1pt} AG)$ является гораздо менее жестким, нежели требование цикличности матрицы A. В-четвертых, хорошо выбранная матрица G должна обеспечивать компромисс между ее размерами и трудностью синтеза стабилизирующего управления.

Если число m является относительно небольшим, то шаг 3 не вызывает серьезных вычислительных проблем. Однако часть вычислительных затрат кроется в процессе предварительной подготовки решения алгебраического уравнения Риккати (3.4). На самом деле хорошо известно, что любой метод решения уравнения Риккати затратнее, чем решение линейных уравнений (0.1). Тем не менее, как указано в [1], что вариации LQRES, например, с помощью субоптимальных методов, использованных для решения уравнения (3.4), может быть более привлекательной численной альтернативой, хорошо известной решателям линейных уравнений. Одним из возможных подходов в этом направлении может быть модельное прогностическое управление, где решение алгебраического уравнения Риккати Р заменяется решением P_m после M шагов разностного уравнения Риккати. Предварительные численные эксперименты с методом LQRESD(M) показали оптимистические результаты по сходимости и устойчивости [1].

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

4. Линейное управление итерационными решателями на основе метода многошаговой декомпозиции. Рассмотрим далее метод линейного управления линейных многомерных динамических систем, используемый для синтеза быстрых и сверхбыстрых итерационных решателей СЛАУ (0.1).

Пусть, как и ранее, $\tilde {A} = {{E}_{n}} - A$, $\tilde {G} = - AG$. Обозначим через ${{\tilde {G}}_{ \bot }}$ матрицу, которая представляет собой левый делитель матрицы $\tilde {G}$ и соответственно удовлетворяет условиям [7]

(4.1)

${{\tilde {G}}_{ \bot }}\tilde {G} = {\mathbf{0}},\quad {{\tilde {G}}_{ \bot }}\tilde {G}_{ \bot }^{{\text{ + }}} = {{E}_{m}},\quad {{\tilde {G}}_{ \bot }}\tilde {G}_{ \bot }^{{\text{ + }}} = {{({{\tilde {G}}_{ \bot }}\tilde {G}_{ \bot }^{{\text{ + }}})}^{{\text{т }}}},\quad \tilde {G}_{ \bot }^{{\text{ + }}}{{\tilde {G}}_{ \bot }} = {{(\tilde {G}_{ \bot }^{{\text{ + }}}{{\tilde {G}}_{ \bot }})}^{{\text{т }}}},$

где $\tilde {G}_{ \bot }^{{\text{ + }}}$ – псевдообратная матрица для матрицы ${{\tilde {G}}_{ \bot }}$.

Определим матрицу регулятора $K$ в формуле закона управления (1.8) в следующем виде:

(4.2)

$K = {{\tilde {G}}^{ + }}\tilde {A} - \Lambda {{\tilde {G}}^{ + }}.$

Здесь ${{\tilde {G}}^{ + }}$ – псевдообратная матрица для матрицы $\tilde {G}$, $\Lambda $ – матрица, имеющая заданные (желаемые) собственные значения.

Можно показать, что справедливым является равенство

(4.3)

${\text{eig(}}\tilde {A} + \tilde {G}K{\text{)}} = {\text{eig(}}{{\tilde {G}}_{ \bot }}\tilde {A}\tilde {G}_{ \bot }^{{\text{ + }}}{\text{)}} \cup {\text{eig(}}\Lambda {\text{)}}{\kern 1pt} {\kern 1pt} .$

Таким образом, закон управления (1.8) с матрицей регулятора (4.2) обеспечивает замкнутой управлением динамической системе заданное подмножество собственных значений, совпадающее с множеством собственных значений матрицы $\Lambda $. Другая часть собственных значений, как видно из формулы (4.3), определяется матрицей ${{\tilde {G}}_{ \bot }}\tilde {A}\tilde {G}_{ \bot }^{{\text{ + }}}$.

Недостаток рассматриваемого закона управления состоит в том, что в общем случае нельзя гарантировать устойчивость матрицы ${{\tilde {G}}_{ \bot }}\tilde {A}\tilde {G}_{ \bot }^{{\text{ + }}}$, т.е. множество ${\text{eig(}}{{\tilde {G}}_{ \bot }}\tilde {A}\tilde {G}_{ \bot }^{{\text{ + }}}{\text{)}}$ может содержать неустойчивые собственные значения ${{\lambda }_{i}}$, у которых ${\text{|}}{{\lambda }_{i}}{\text{|}} > 1$. Поэтому введем в рассмотрение многошаговую декомпозицию системы (3.1). Декомпозиция имеет следующий вид:

исходное положение

(4.4)

${{\tilde {A}}_{0}} = \tilde {A},\quad {{\tilde {G}}_{0}} = \tilde {G},$

первый шаг

(4.5)

${{\tilde {A}}_{1}} = {{\tilde {G}}_{ \bot }}\tilde {A}\tilde {G}_{ \bot }^{{\text{ + }}},\quad {{\tilde {G}}_{1}} = {{\tilde {G}}_{ \bot }}\tilde {A}\tilde {G},$

k-й (промежуточный) шаг

(4.6)

${{\tilde {A}}_{k}} = {{\tilde {G}}_{{ \bot k - 1}}}{{\tilde {A}}_{{k - 1}}}\tilde {G}_{{ \bot k - 1}}^{{\text{ + }}},\quad {{\tilde {G}}_{k}} = {{\tilde {G}}_{{ \bot k - 1}}}{{\tilde {A}}_{{k - 1}}}{{\tilde {G}}_{{k - 1}}},$

S-й (конечный) шаг

(4.7)

${{\tilde {A}}_{S}} = {{\tilde {G}}_{{ \bot S - 1}}}{{\tilde {A}}_{{S - 1}}}\tilde {G}_{{ \bot S - 1}}^{{\text{ + }}},\quad {{\tilde {G}}_{S}} = {{\tilde {G}}_{{ \bot S - 1}}}{{\tilde {A}}_{{S - 1}}}{{\tilde {G}}_{{S - 1}}}.$

Здесь $S = {\text{floor(}}n{\text{/}}m{\text{)}} - 1$; floor – операция округления числа $n{\text{/}}m$ в сторону ближайшего целого в меньшую сторону, например, ${\text{floor}}(0.4) = 0$, ${\text{floor}}(2.9) = 2$ и т.д.

Если также пошагово определить матрицы

(4.8)

${{K}_{0}} = {{\Lambda }_{0}}\tilde {G}_{0}^{ - } - \tilde {G}_{0}^{ - }{{\tilde {A}}_{0}}{\text{,}}\quad \tilde {G}_{0}^{ - } = \tilde {G}_{0}^{ + } - {{K}_{1}}{{\tilde {G}}_{{ \bot 0}}},$

(4.9)

${{K}_{1}} = {{\Lambda }_{1}}\tilde {G}_{1}^{ - } - \tilde {G}_{1}^{ - }{{\tilde {A}}_{1}}{\text{,}}\quad \tilde {G}_{1}^{ - } = \tilde {G}_{1}^{ + } - {{K}_{2}}{{\tilde {G}}_{{ \bot 1}}},\; \ldots ,$

(4.10)

${{K}_{k}} = {{\Lambda }_{k}}\tilde {G}_{k}^{ - } - \tilde {G}_{k}^{ - }{{\tilde {A}}_{k}}{\text{,}}\quad \tilde {G}_{k}^{ - } = \tilde {G}_{k}^{ + } - {{K}_{{k + 1}}}{{\tilde {G}}_{{ \bot k}}},\; \ldots ,$

(4.11)

${{K}_{S}} = {{\Lambda }_{S}}\tilde {G}_{S}^{ - } - \tilde {G}_{S}^{ + }{{\tilde {A}}_{S}},$

то в результате замыкания обратной связи с законом управления (1.8), где матрица $K$ имеет вид

(4.12)

$K = \Lambda ({{\tilde {G}}^{ + }} - {{K}_{1}}{{\tilde {G}}_{ \bot }}) - ({{\tilde {G}}^{ + }} - {{K}_{1}}{{\tilde {G}}_{ \bot }})\tilde {A},$

обеспечивается следующее равенство множеств собственных значений [7]:

(4.13)

${\text{eig(}}\tilde {A} - \tilde {G}K{\text{)}} = \bigcup\limits_{i = 1}^{S + 1} {{\text{eig}}({{\Lambda }_{{i - 1}}})} .$

Здесь ${{\Lambda }_{{i - 1}}}$, $i - 1 = 0,1, \ldots ,S$, – матрицы, имеющие заданные собственные значения на каждом уровне декомпозиции, причем

$\bigcup\limits_{i = 1}^{S + 1} {{\text{eig}}({{\Lambda }_{{i - 1}}}) = {\text{eig}}(\Lambda ).} $

Таким образом, многошаговая декомпозиция системы (4.4)–(4.7) и такая же по количеству шагов процедура определения матрицы K (4.12) на основе выражений (4.8)–(4.11) обеспечивают заданное размещение собственных значений $\Lambda = \{ {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\lambda } }_{1}},{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\lambda } }_{2}}, \ldots {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\lambda } }_{n}}\} $, как это указано в формуле (4.13).

Рассмотрим задачу синтеза закона управления (1.8), обеспечивающего конечную длительность переходного процесса дискретной MIMO-системы. В этом случае матрица $\tilde {A} - \tilde {G}K$ в качестве собственных значений должна иметь только нули [26], т.е.

(4.14)

${\text{eig}}(\tilde {A} + \tilde {G}K) = \{ 0,0, \ldots ,0\} .$

Это требование означает, что в качестве матриц ${{\Lambda }_{i}}$ могут быть взяты любые нильпотентные матрицы с индексом нильпотентности не больше $r$ [4, 7]. Тогда сходимость синтезированных решателей, рассматриваемых как дискретные системы, будет наступать не более чем за n шагов. Такие решатели будем называть быстрыми.

Если в качестве матриц ${{\Lambda }_{i}}$ выбрать нулевые матрицы, то этому случаю можно условно сопоставить дискретную систему ${{x}_{{k + 1}}} = b$, сходимость которой, очевидно, наступает на втором шаге, т.е. ${{x}_{0}} = 0$, ${{x}_{1}} = b$. Подобные решатели можно охарактеризовать как сверхбыстрые.

При нулевых матрицах соответствующего размера ${{\Lambda }_{i}}$ заключительная часть алгоритма (4.8)–(4.11) примет вид

(4.15)

${{K}_{0}} = - (\tilde {G}_{0}^{ + } - {{K}_{1}}{{\tilde {G}}_{{ \bot 0}}}){{\tilde {A}}_{0}}{\text{,}}$

(4.16)

${{K}_{1}} = - (\tilde {G}_{1}^{ + } - {{K}_{2}}{{\tilde {G}}_{{ \bot 1}}}){{\tilde {A}}_{1}}{\text{,}} \ldots {\text{,}}$

(4.17)

${{K}_{k}} = - (\tilde {G}_{k}^{ + } - {{K}_{{k + 1}}}{{\tilde {G}}_{{ \bot k}}}){{\tilde {A}}_{k}}{\text{,}} \ldots {\text{,}}$

(4.18)

${{K}_{S}} = - \tilde {G}_{S}^{ + }{{\tilde {A}}_{S}}.$

Как следует из (4.15)–(4.18), при количестве шагов декомпозиции $S = 1$ формула матрицы K в законе управления (1.8), обеспечивающего не только конечную длительность переходного процесса, но и сверхбыструю (в течение первых итераций независимо от размерности дискретной MIMO-системы) сходимость, принимает следующий вид:

(4.19)

$K = - ({{\tilde {G}}^{ + }} + \tilde {G}_{1}^{ + }{{\tilde {A}}_{1}}{{\tilde {G}}_{ \bot }})\tilde {A}.$

При $S = 2$ и $S = 3$ формулы матрицы K несколько усложняются:

(4.20)

$K = - ({{\tilde {G}}^{ + }} + (\tilde {G}_{1}^{ + } + \tilde {G}_{2}^{ + }{{\tilde {A}}_{2}}{{\tilde {G}}_{{ \bot 1}}}){{\tilde {A}}_{1}}{{\tilde {G}}_{ \bot }})\tilde {A},$

(4.21)

$K = - ({{\tilde {G}}^{ + }} + (\tilde {G}_{1}^{ + } + (\tilde {G}_{2}^{ + } + \tilde {G}_{3}^{ + }{{\tilde {A}}_{3}}{{\tilde {G}}_{{ \bot 2}}}){{\tilde {A}}_{2}}{{\tilde {G}}_{{ \bot 1}}}){{\tilde {A}}_{1}}{{\tilde {G}}_{ \bot }})\tilde {A}.$

Как следует из анализа конечных формул (4.19)–(4.21), синтезированные на основе многошаговой декомпозиции матрицы K содержат вложения, вид которых определяется как собственно матрицей A системы (0.1), так и “произвольно” выбираемой матрицей G. При этом отсутствие ограничений размеров матрицы G предоставляет широкие возможности по выбору матриц K. Так, если заранее определить количество столбцов матрицы G не менее чем n/2, то закон управления с матрицей K (4.19) обеспечивает быструю сходимость итерационного процесса следующего вида:

(4.22)

${{x}_{{k + 1}}} = \tilde {A}{{x}_{k}} - G({{\tilde {G}}^{ + }} + \tilde {G}_{1}^{ + }{{\tilde {A}}_{1}}{{\tilde {G}}_{ \bot }})\tilde {A}(b - A{{x}_{k}}) + b.$

5. Примеры синтеза сверхбыстрых итерационных решателей. Рассмотрим сначала СЛАУ (0.1) в общем виде с нециклической и не удовлетворяющей условиям нормальности матрицей:

(5.1)

$A = \left[ {\begin{array}{*{20}{c}} r&w&\vline & 0&0 \\ { - w}&r&\vline & 0&0 \\ \hline 0&0&\vline & \lambda &1 \\ 0&0&\vline & 0&\lambda \end{array}} \right] \in {{\mathbb{R}}^{{4 \times 4}}},$

где r, w, λ – некоторые заданные числа. На нецикличность матрицы (5.1) указывают диагональные блоки 2 × 2, где в первом в явном виде записана ячейка Жордана, а во втором – верхнетреугольная матрица. Ненормальность данной матрицы проверяется непосредственно: ${{A}^{{\text{т }}}}A \ne A{{A}^{{\text{т }}}}$.

Пусть вектор b также задан в общем виде

(5.2)

$b = {{[\begin{array}{*{20}{c}} {b{}_{1}}&{b{}_{2}}&{b{}_{3}}&{b{}_{4}} \end{array}]}^{{\text{т }}}}.$

Зададим матрицу G, например, в следующем виде:

(5.3)

$G = {{\left[ {\begin{array}{*{20}{c}} 1&0&0&0 \\ 0&0&0&1 \end{array}} \right]}^{{\text{т }}}} \in {{\mathbb{R}}^{{4 \times 2}}}.$

Выполним синтез итерационного решателя согласно алгоритму (4.4)–(4.11), (4.15)–(4.18). Как видно из анализа размеров матриц (5.1) и (5.3), в данном случае S = floor(n/m) – 1 = ${\text{floor}}(4{\text{/}}2) - 1 = 2 - 1 = 1$ и, следовательно, выполнив двухшаговую декомпозицию, включающую нулевой и первый шаг, можно воспользоваться формулой для матрицы K (4.19).

Итак, в соответствии с выражениями (4.7) и (4.8) вычислим матрицы ${{\tilde {A}}_{0}}$, ${{\tilde {G}}_{0}}$, ${{\tilde {A}}_{1}}$, ${{\tilde {G}}_{1}}$. Имеем

(5.4)

${{\tilde {A}}_{0}} = \left[ {\begin{array}{*{20}{c}} {1 - r}&{ - w}&\vline & 0&0 \\ w&{1 - r}&\vline & 0&0 \\ \hline 0&0&\vline & {1 - \lambda }&{ - 1} \\ 0&0&\vline & 0&{1 - \lambda } \end{array}} \right],\quad {{\tilde {G}}_{0}} = \left[ {\begin{array}{*{20}{c}} { - r}&0 \\ w&0 \\ 0&{ - 1} \\ 0&{ - \lambda } \end{array}} \right],$

(5.5)

${{\tilde {A}}_{1}} = \left[ {\begin{array}{*{20}{c}} {1 - r}&0 \\ 0&{\frac{{ - {{\lambda }^{3}} + {{\lambda }^{2}} + 1}}{{{{\lambda }^{2}} + 1}}} \end{array}} \right],\quad {{\tilde {G}}_{1}} = \left[ {\begin{array}{*{20}{c}} { - \frac{{w({{r}^{2}} + {{w}^{2}})}}{r}}&0 \\ 0&{ - {{\lambda }^{2}}} \end{array}} \right].$

Отметим, что для вычисления матриц (5.5) использовались следующие матрицы:

(5.6)

${{\tilde {G}}_{ \bot }} = \left[ {\begin{array}{*{20}{c}} {\frac{w}{r}}&1&0&0 \\ 0&0&{ - \lambda }&1 \end{array}} \right]{\text{,}}\quad \tilde {G}_{ \bot }^{{\text{ + }}} = \left[ {\begin{array}{*{20}{c}} {\frac{{rw}}{{{{r}^{2}} + {{w}^{2}}}}}&{\frac{{{{r}^{2}}}}{{{{r}^{2}} + {{w}^{2}}}}}&0&0 \\ 0&0&{ - \frac{\lambda }{{{{\lambda }^{2}} + 1}}}&{\frac{1}{{{{\lambda }^{2}} + 1}}} \end{array}} \right].$

Как следует из (4.19), для определения матрицы K также необходима обратная матрица $\tilde {G}_{1}^{ + }$, которая в соответствии с (5.5) равна

(5.7)

$\tilde {G}_{1}^{ + } = \left[ {\begin{array}{*{20}{c}} { - \frac{r}{{w({{r}^{2}} + {{w}^{2}})}}}&0 \\ 0&{ - \frac{1}{{{{\lambda }^{2}}}}} \end{array}} \right].$

Подставляя матрицы $\tilde {A} = {{\tilde {A}}_{0}}$ (5.4), ${{\tilde {A}}_{1}}$ (5.5), ${{\tilde {G}}_{ \bot }}$ (5.6) и $\tilde {G}_{1}^{ + }$ (5.7) в формулу K (4.19), получим

$K = \left[ {\begin{array}{*{20}{c}} { - \frac{{{{r}^{2}} + {{w}^{2}} - 1}}{{{{r}^{2}} + {{w}^{2}}}}}&{\frac{{{{r}^{3}} - 2{{r}^{2}} + r{{w}^{2}} + r - 2{{w}^{2}}}}{{{{r}^{2}} + {{w}^{2}}}}}&0&0 \\ 0&0&{ - \frac{{{{{(\lambda - 1)}}^{2}}}}{\lambda }}&{ - \frac{{{{\lambda }^{2}} - 1}}{{{{\lambda }^{2}}}}} \end{array}} \right].$

При этом

$\tilde {A} + \tilde {G}K = \left[ {\begin{array}{*{20}{c}} {\frac{{{{r}^{2}} - r + {{w}^{2}}}}{{{{r}^{2}} + {{w}^{2}}}}}&{ - \frac{{{{{({{r}^{2}} - r + {{w}^{2}})}}^{2}}}}{{w({{r}^{2}} + {{w}^{2}})}}}&\vline & 0&0 \\ {\frac{w}{{{{r}^{2}} + {{w}^{2}}}}}&{ - \frac{{{{r}^{2}} - r + {{w}^{2}}}}{{{{r}^{2}} + {{w}^{2}}}}}&\vline & 0&0 \\ \hline 0&0&\vline & { - \frac{{\lambda - 1}}{\lambda }}&{ - \frac{1}{{{{\lambda }^{2}}}}} \\ 0&0&\vline & {{{{(\lambda - 1)}}^{2}}}&{\frac{{\lambda - 1}}{\lambda }} \end{array}} \right],$

${\text{eig(}}\tilde {A} + \tilde {G}K{\text{)}} = \{ 0,0,0,0\} ,$

что и требовалось.

Запуская далее синтезированный итерационный решатель (4.22) с указанными ранее параметрами для нулевых начальных условий, получим следующую последовательность невязок:

${{e}_{0}} = {{[\begin{array}{*{20}{c}} {b{}_{1}}&{b{}_{2}}&{b{}_{3}}&{b{}_{4}} \end{array}]}^{{\text{т }}}},$

${{e}_{1}} = \left[ {\begin{array}{*{20}{c}} { - \frac{{({{r}^{2}} - r + {{w}^{2}})( - {{b}_{2}}{{r}^{2}} + {{b}_{2}}r - {{b}_{2}}{{w}^{2}} + {{b}_{1}}w)}}{{w({{r}^{2}} + {{w}^{2}})}}} \\ { - \frac{{ - {{b}_{2}}{{r}^{2}} + {{b}_{2}}r - {{b}_{2}}{{w}^{2}} + {{b}_{1}}w}}{{{{r}^{2}} + {{w}^{2}}}}} \\ {\frac{{{{b}_{3}}{{\lambda }^{2}} - {{b}_{3}}\lambda + {{b}_{4}}}}{{{{\lambda }^{2}}}}} \\ { - \frac{{(\lambda - 1)({{b}_{3}}{{\lambda }^{2}} - {{b}_{3}}\lambda + {{b}_{4}})}}{\lambda }} \end{array}} \right],$

${{e}_{2}} = {{[\begin{array}{*{20}{c}} 0&0&0&0 \end{array}]}^{{\text{T}}}}, \ldots ,{{e}_{k}} = {{[\begin{array}{*{20}{c}} 0&0&0&0 \end{array}]}^{{\text{T}}}},\quad k > 2.$

Как видно, итерационный процесс сошелся к точному решению уже на второй итерации.

Такая же сходимость характерна, например, для следующей трехдиагональной матрицы A, матрицы G и вектора b:

(5.8)

$A = \left[ {\begin{array}{*{20}{c}} {{{a}_{{11}}}}&{{{a}_{{12}}}}&0&0 \\ {{{a}_{{21}}}}&{{{a}_{{22}}}}&{{{a}_{{23}}}}&0 \\ 0&{{{a}_{{32}}}}&{{{a}_{{33}}}}&{{{a}_{{34}}}} \\ 0&0&{{{a}_{{43}}}}&{{{a}_{{44}}}} \end{array}} \right] \in {{\mathbb{R}}^{{4 \times 4}}},\quad G = \left[ {\begin{array}{*{20}{c}} 0&0 \\ 1&0 \\ 0&0 \\ 0&1 \end{array}} \right] \in {{\mathbb{R}}^{{4 \times 2}}},\quad b = \left[ {\begin{array}{*{20}{c}} {{{b}_{1}}} \\ {{{b}_{2}}} \\ {{{b}_{3}}} \\ {{{b}_{4}}} \end{array}} \right] \in {{\mathbb{R}}^{4}}.$

Продемонстрируем это. Как и прежде, согласно выражениям (4.7) и (4.8), определим матрицы ${{\tilde {A}}_{0}}$, ${{\tilde {G}}_{0}}$, ${{\tilde {A}}_{1}}$, ${{\tilde {G}}_{1}}$:

(5.9)

${{\tilde {A}}_{0}} = \left[ {\begin{array}{*{20}{c}} {1 - {{a}_{{11}}}}&{ - {{a}_{{12}}}}&0&0 \\ { - {{a}_{{21}}}}&{1 - {{a}_{{22}}}}&{ - {{a}_{{23}}}}&0 \\ 0&{ - {{a}_{{32}}}}&{1 - {{a}_{{33}}}}&{ - {{a}_{{34}}}} \\ 0&0&{ - {{a}_{{43}}}}&{1 - {{a}_{{44}}}} \end{array}} \right],\quad {{\tilde {G}}_{0}} = \left[ {\begin{array}{*{20}{c}} { - {{a}_{{12}}}}&0 \\ { - {{a}_{{22}}}}&0 \\ { - {{a}_{{32}}}}&{ - {{a}_{{34}}}} \\ 0&{ - {{a}_{{44}}}} \end{array}} \right],$

(5.10)

${{\tilde {A}}_{1}} = \frac{1}{\Delta }\left[ {\begin{array}{*{20}{c}} {{{{\tilde {a}}}_{{11}}}}&{{{{\tilde {a}}}_{{12}}}} \\ {{{{\tilde {a}}}_{{21}}}}&{{{{\tilde {a}}}_{{22}}}} \end{array}} \right],\quad {{\tilde {G}}_{1}} = \left[ {\begin{array}{*{20}{c}} {{{a}_{{12}}}{{a}_{{21}}} - {{a}_{{11}}}{{a}_{{22}}} + {{a}_{{23}}}{{a}_{{32}}}}&{{{a}_{{23}}}{{a}_{{34}}}} \\ {\frac{{{{a}_{{32}}}({{a}_{{11}}}{{a}_{{44}}} - {{a}_{{33}}}{{a}_{{44}}} + {{a}_{{34}}}{{a}_{{43}}})}}{{{{a}_{{34}}}}}}&{{{a}_{{34}}}{{a}_{{43}}} - {{a}_{{33}}}{{a}_{{44}}}} \end{array}} \right],$

где для компактности записи ведены обозначения:

При вычислении матриц (5.10) применялись следующие матрицы:

(5.11)

${{\tilde {G}}_{ \bot }} = \left[ {\begin{array}{*{20}{c}} { - \frac{{{{a}_{{22}}}}}{{{{a}_{{12}}}}}}&1&0&0 \\ {\frac{{{{a}_{{32}}}{{a}_{{44}}}}}{{{{a}_{{12}}}{{a}_{{34}}}}}}&0&{ - \frac{{{{a}_{{44}}}}}{{{{a}_{{34}}}}}}&1 \end{array}} \right],\quad G_{ \bot }^{{\text{ + }}} = \frac{1}{\Delta }\left[ {\begin{array}{*{20}{c}} { - {{a}_{{12}}}{{a}_{{22}}}({{a}_{{34}}}^{2} + {{a}_{{44}}}^{2})}&{{{a}_{{12}}}{{a}_{{32}}}{{a}_{{34}}}{{a}_{{44}}}} \\ {{{a}_{{12}}}^{2}{{a}_{{34}}}^{2} + {{a}_{{12}}}^{2}{{a}_{{44}}}^{2} + {{a}_{{32}}}^{2}{{a}_{{44}}}^{2}}&{{{a}_{{22}}}{{a}_{{32}}}{{a}_{{34}}}{{a}_{{44}}}} \\ { - {{a}_{{22}}}{{a}_{{32}}}{{a}_{{44}}}^{2}}&{ - {{a}_{{34}}}{{a}_{{44}}}({{a}_{{12}}}^{2} + {{a}_{{22}}}^{2})} \\ {{{a}_{{22}}}{{a}_{{32}}}{{a}_{{34}}}{{a}_{{44}}}}&{{{a}_{{34}}}^{2}({{a}_{{12}}}^{2} + {{a}_{{22}}}^{2})} \end{array}} \right].$

При этом псевдообратные матрицы $\tilde {G}_{{}}^{ + }$, $\tilde {G}_{1}^{ + }$ равны

(5.12)

$\tilde {G}_{{}}^{ + } = \frac{1}{\Delta }{{\left[ {\begin{array}{*{20}{c}} { - {{a}_{{12}}}({{a}_{{34}}}^{2} + {{a}_{{44}}}^{2})}&{{{a}_{{12}}}{{a}_{{32}}}{{a}_{{34}}}} \\ { - {{a}_{{22}}}({{a}_{{34}}}^{2} + a{{{_{{44}}^{{}}}}^{2}})}&{{{a}_{{22}}}{{a}_{{32}}}{{a}_{{34}}}} \\ { - {{a}_{{32}}}{{a}_{{44}}}^{2}}&{ - {{a}_{{34}}}({{a}_{{12}}}^{2} + {{a}_{{22}}}^{2})} \\ {{{a}_{{32}}}{{a}_{{34}}}{{a}_{{44}}}}&{ - {{a}_{{44}}}({{a}_{{12}}}^{2} + {{a}_{{22}}}^{2} + {{a}_{{32}}}^{2})} \end{array}} \right]}^{{\text{т }}}},$

(5.13)

$\tilde {G}_{1}^{ + } = \frac{1}{{{{\Delta }_{1}}}}\left[ {\begin{array}{*{20}{c}} {{{a}_{{33}}}{{a}_{{44}}} - {{a}_{{34}}}{{a}_{{43}}}}&{{{a}_{{23}}}{{a}_{{34}}}} \\ {{{a}_{{32}}}({{a}_{{11}}}{{a}_{{44}}} - {{a}_{{33}}}{{a}_{{44}}} + {{a}_{{34}}}{{a}_{{43}}})}&{ - ({{a}_{{12}}}{{a}_{{21}}} - {{a}_{{11}}}{{a}_{{22}}} + {{a}_{{23}}}{{a}_{{32}}})} \end{array}} \right],$

где ${{\Delta }_{1}} = {{a}_{{11}}}{{a}_{{22}}}{{a}_{{34}}}{{a}_{{43}}} - {{a}_{{11}}}{{a}_{{22}}}{{a}_{{33}}}a_{{44}}^{{}} + {{a}_{{11}}}{{a}_{{23}}}{{a}_{{32}}}{{a}_{{44}}} + {{a}_{{12}}}{{a}_{{21}}}{{a}_{{33}}}{{a}_{{44}}} - {{a}_{{12}}}{{a}_{{21}}}{{a}_{{34}}}{{a}_{{43}}}$.

Используя далее матрицы $\tilde {A} = {{\tilde {A}}_{0}}$ (5.9), ${{\tilde {A}}_{1}}$ (5.10), ${{\tilde {G}}_{ \bot }}$ (5.11), $\tilde {G}_{{}}^{ + }$ (5.12), $\tilde {G}_{1}^{ + }$ (5.13) и формулу для K (4.19), получим

$K = \frac{1}{{{{\Delta }_{K}}}}\left[ {\begin{array}{*{20}{c}} {{{k}_{{11}}}}&{{{k}_{{12}}}}&{{{k}_{{13}}}}&{{{k}_{{14}}}} \\ {{{k}_{{21}}}}&{{{k}_{{22}}}}&{{{k}_{{23}}}}&{{{k}_{{24}}}} \end{array}} \right].$

Здесь введены обозначения:

${{\Delta }_{K}} = {{a}_{{12}}}({{a}_{{11}}}{{a}_{{22}}}{{a}_{{34}}}{{a}_{{43}}} - {{a}_{{11}}}{{a}_{{22}}}{{a}_{{33}}}{{a}_{{44}}} + {{a}_{{11}}}{{a}_{{23}}}{{a}_{{32}}}{{a}_{{44}}} + {{a}_{{12}}}{{a}_{{21}}}{{a}_{{33}}}{{a}_{{44}}} - {{a}_{{12}}}{{a}_{{21}}}{{a}_{{34}}}{{a}_{{43}}}),$

$\begin{gathered} {{k}_{{11}}} = - ({{a}_{{22}}}{{a}_{{34}}}{{a}_{{43}}} - {{a}_{{22}}}{{a}_{{33}}}{{a}_{{44}}} + {{a}_{{23}}}{{a}_{{32}}}{{a}_{{44}}} - {{a}_{{11}}}^{2}{{a}_{{22}}}{{a}_{{33}}}{{a}_{{44}}} + {{a}_{{11}}}^{2}{{a}_{{22}}}{{a}_{{34}}}{{a}_{{43}}} + \\ \, + {{a}_{{11}}}^{2}{{a}_{{23}}}{{a}_{{32}}}{{a}_{{44}}} + 2{{a}_{{11}}}{{a}_{{22}}}{{a}_{{33}}}{{a}_{{44}}} - 2{{a}_{{11}}}{{a}_{{22}}}{{a}_{{34}}}{{a}_{{43}}} - 2{{a}_{{11}}}{{a}_{{23}}}{{a}_{{32}}}{{a}_{{44}}} - 2{{a}_{{12}}}{{a}_{{21}}}{{a}_{{33}}}{{a}_{{44}}} + \\ \, + 2{{a}_{{12}}}{{a}_{{21}}}{{a}_{{34}}}{{a}_{{43}}} + {{a}_{{11}}}{{a}_{{12}}}{{a}_{{21}}}{{a}_{{33}}}{{a}_{{44}}} - {{a}_{{11}}}{{a}_{{12}}}{{a}_{{21}}}{{a}_{{34}}}{{a}_{{43}}}), \\ \end{gathered} $

$\begin{gathered} {{k}_{{12}}} = - ({{a}_{{33}}}{{a}_{{44}}} - {{a}_{{34}}}{{a}_{{43}}} - {{a}_{{11}}}{{a}_{{22}}}{{a}_{{33}}}{{a}_{{44}}} + {{a}_{{11}}}{{a}_{{22}}}{{a}_{{34}}}{{a}_{{43}}} + {{a}_{{11}}}{{a}_{{23}}}{{a}_{{32}}}{{a}_{{44}}} + \\ \, + {{a}_{{12}}}{{a}_{{21}}}{{a}_{{33}}}{{a}_{{44}}} - a{}_{{12}}{{a}_{{21}}}{{a}_{{34}}}{{a}_{{43}}})a_{{21}}^{{ - 1}}, \\ \end{gathered} $

${{k}_{{13}}} = {{a}_{{23}}}{{a}_{{44}}},\quad {{k}_{{14}}} = - {{a}_{{23}}}{{a}_{{34}}},$

$\begin{gathered} {{k}_{{21}}} = {{a}_{{32}}}({{a}_{{12}}}{{a}_{{21}}}{{a}_{{44}}} - {{a}_{{22}}}{{a}_{{33}}}{{a}_{{44}}} + {{a}_{{22}}}{{a}_{{34}}}{{a}_{{43}}} + a{}_{{23}}{{a}_{{32}}}{{a}_{{44}}} - {{a}_{{11}}}^{2}{{a}_{{22}}}{{a}_{{33}}}{{a}_{{44}}} + \\ \, + {{a}_{{11}}}^{2}{{a}_{{22}}}{{a}_{{34}}}{{a}_{{43}}} + {{a}_{{11}}}^{2}{{a}_{{23}}}{{a}_{{32}}}{{a}_{{44}}} + 2{{a}_{{11}}}{{a}_{{22}}}{{a}_{{33}}}{{a}_{{44}}} - 2{{a}_{{11}}}{{a}_{{22}}}{{a}_{{34}}}{{a}_{{43}}} - 2{{a}_{{11}}}{{a}_{{23}}}{{a}_{{32}}}{{a}_{{44}}} - \\ \, - 2{{a}_{{12}}}{{a}_{{21}}}{{a}_{{33}}}{{a}_{{44}}} + 2{{a}_{{12}}}{{a}_{{21}}}{{a}_{{34}}}{{a}_{{43}}} + {{a}_{{11}}}{{a}_{{12}}}{{a}_{{21}}}{{a}_{{33}}}{{a}_{{44}}} - {{a}_{{11}}}{{a}_{{12}}}{{a}_{{21}}}{{a}_{{34}}}{{a}_{{43}}}){{({{a}_{{12}}}{{a}_{{34}}})}^{{ - 1}}}, \\ \end{gathered} $

${{k}_{{22}}} = - {{a}_{{32}}}({{a}_{{11}}}{{a}_{{44}}} - {{a}_{{33}}}{{a}_{{44}}} + {{a}_{{34}}}{{a}_{{43}}})a_{{34}}^{{ - 1}},$

$\begin{gathered} {{k}_{{23}}} = - ({{a}_{{12}}}{{a}_{{21}}}{{a}_{{44}}} - {{a}_{{11}}}{{a}_{{22}}}{{a}_{{44}}} + {{a}_{{23}}}{{a}_{{32}}}{{a}_{{44}}} - {{a}_{{11}}}{{a}_{{22}}}{{a}_{{33}}}^{2}{{a}_{{44}}} + {{a}_{{12}}}{{a}_{{21}}}{{a}_{{33}}}^{2}{{a}_{{44}}} + \\ \, + 2{{a}_{{11}}}{{a}_{{22}}}{{a}_{{33}}}{{a}_{{44}}} - 2{{a}_{{11}}}{{a}_{{22}}}a{}_{{34}}{{a}_{{43}}} - 2{{a}_{{11}}}{{a}_{{23}}}{{a}_{{32}}}{{a}_{{44}}} - 2{{a}_{{12}}}{{a}_{{21}}}{{a}_{{33}}}{{a}_{{44}}} + 2{{a}_{{12}}}{{a}_{{21}}}{{a}_{{34}}}{{a}_{{43}}} + \\ \, + {{a}_{{11}}}{{a}_{{22}}}{{a}_{{33}}}{{a}_{{34}}}{{a}_{{43}}} + {{a}_{{11}}}{{a}_{{23}}}{{a}_{{32}}}{{a}_{{33}}}{{a}_{{44}}} - {{a}_{{12}}}{{a}_{{21}}}{{a}_{{33}}}{{a}_{{34}}}{{a}_{{43}}})a_{{34}}^{{ - 1}}, \\ \end{gathered} $

$\begin{gathered} {{k}_{{24}}} = {{a}_{{11}}}{{a}_{{22}}} - {{a}_{{12}}}{{a}_{{21}}} - {{a}_{{23}}}{{a}_{{32}}} - {{a}_{{11}}}{{a}_{{22}}}{{a}_{{33}}}{{a}_{{44}}} + {{a}_{{11}}}{{a}_{{22}}}{{a}_{{34}}}{{a}_{{43}}} + {{a}_{{11}}}{{a}_{{23}}}{{a}_{{32}}}{{a}_{{44}}} + {\text{ }} \\ \, + {{a}_{{12}}}{{a}_{{21}}}{{a}_{{33}}}{{a}_{{44}}} - {{a}_{{12}}}{{a}_{{21}}}{{a}_{{34}}}{{a}_{{43}}}. \\ \end{gathered} $

Как и в предыдущем случае, при начальной невязке ${{e}_{0}} = {{[\begin{array}{*{20}{c}} {b{}_{1}}&{b{}_{2}}&{b{}_{3}}&{b{}_{4}} \end{array}]}^{{\text{т }}}}$ вычислительный процесс определяет точное решение на второй итерации.

Исследование синтезированных сверхбыстрых решателей линейных уравнений с матрицами и векторами, элементы которых подчиняются нормальному закону распределения, а собственные значения матриц A – круговому закону Гирко [27]:

$\begin{gathered} A = {\text{randn}}(100 \times 100),\quad G = {\text{randn}}(100 \times 50),\quad b = {\text{randn}}(100 \times 1), \\ A = {\text{randn}}(500 \times 500),\quad G = {\text{randn}}(500 \times 250),\quad b = {\text{randn}}(500 \times 1), \\ A = {\text{randn}}(1000 \times 1000),\quad G = {\text{randn}}(1000 \times 500),\quad b = {\text{randn}}(1000 \times 1), \\ A = {\text{randn}}(5000 \times 5000),\quad G = {\text{randn}}(5000 \times 2500),\quad b = {\text{randn}}(5000 \times 1) \\ \end{gathered} $

показали, что во всех случаях сходимость итерационного процесса наступала на третьей или четвертой итерации независимо от размерности решаемой задачи. При этом для каждой размерности задачи было проведено не менее n+1 испытаний.

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

Сложность реализации предложенного итерационного алгоритма могла быть обусловлена отсутствием экономных процедур вычисления псевдообратных матриц неполного ранга, однако в соответствии с работами авторов [28, 29] она может быть успешно преодолена.

Применение предлагаемого метода для решения практических задач, в том числе для реальной электроэнергетической системы, содержащей более 100 электрических станций [30–32], из-за ограничения объема настоящей статьи будет рассмотрено авторами в будущих публикациях.

Список литературы

Helmke U., Jordan J. Control and Stabilization of Linear Equation Solvers // Lecture Notes in Control and Information Sciences. Berlin: Springer-Verlag, 2010. P. 73–82.
Helmke U., Jordan J. Optimal Control of Iterative Solution Methods for Linear Systems of Equations // Proc. Appl. Math. Mech. 2005. V. 5(1). P. 163–164.
Helmke U., Jordan J. Mathematical Systems Theory in Biology, Communications, Computations and Finance // IMA Volumes in Mathematics and its Applications. 2002. V. 134. P. 223–236.
Bernstein D.S. Matrix Mathematics. Princeton: Princeton Univ. Press, 2009.
Голуб Дж., Ван Лоун Ч. Матричные вычисления. М.: Мир, 1999.
Embree M. The Tortoise and the Hare Restart GMRES // SIAM Review. 2003. V. 45. P. 259–266.
Зубов Н.Е., Микрин Е.А., Рябченко В.Н. Матричные методы в теории и практике систем автоматического управления летательными аппаратами. М.: Изд-во МГТУ им. Н.Э. Баумана, 2016.
Bellman R.E. Introduction to the Mathematical Theory of Control Processes // Nonlinear Processes. V. II. N. Y.: Acad. Press, 1971.
Gustafsson K., Lundh M., Soderlind G.S. A PI Stepsize Control for the Numerical Solution of Ordinary Differential Equations // BIT. 1988. V. 28. P. 270–287.
Gustafsson K. Control Theoretic Techniques for Stepsize Selection in Explicit Runge-Kutta Methods // ACM Trans. Math. Softw. 1991. V. 17(4). P. 533–554.
Helmke U., Fuhrmann P.A. Controllability of Matrix Eigenvalue Algorithms: the Inverse Power Method // Systems & Control Letters. 2000. V. 41. P. 57–66.
Helmke U., Jordan J., Lanzon A. A Control Theory Approach to Linear Equation Solvers // Proc. of 17th Intern. Symp. on Mathematical Theory of Network and Systems (MTNS), Kyoto, Japan, 2006.
Helmke U., Wirth F. On Controllability of the Real Shifted Inverse Power Iteration // Systems & Control Letters. 2001. V. 43. P. 9–23.
Jordan J. Reachable Sets of Numerical Iteration Schemes. Ph.D. Thesis. Wurzburg: University of Wurzburg, 2008.
Bhaya A., Kaszkurewicz E. Control Perspectives on Numerical Algorithms and Matrix Problems // Advances in Design and Control. 2006. V. 10.
Kashima K., Yamamoto Y. System Theory for Numerical Analysis // Automatic. 2007. V. 43(7). P. 1156–1164.
Wakasa Y., Yamamoto Y. An Iterative Method for Solving Linear Equations by Using Robust Methods // Proc. of SICE first Annual Conf. on Control Systems. Japan, 2001. P. 451–454.
Calvetti D., Reichel L. An Adaptive Richardson Iteration Method for Indefinite Linear Systems // Numerical Algorithms. 1996. V. 12. P. 125–149.
Golub G.H., Overton M.L. The Convergence of Inexact Chebyshev and Richardson Iterative Methods for Solving Linear Systems // Numer. Math. 1988. V. 53. P. 571–593.
Opfer G., Schober G. Richardson’s Iteration for Nonsymmetric Matrices // Linear Algebra and its Appl. 1984. V. 58. P. 343–361.
Smolarski D.C., Saylor P.E. An Optimum Iterative Method for Solving any Linear Systems with a Square Matrix // BIT. V. 28. P. 163–178.
Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989.
Eiermann M., Ernst O.G., Schneider O. Analysis of Acceleration Strategies for Restarted Minimal Residual Methods // J. Comp. Appl. Math. 2000. V. 123. P. 261–292.
Joubert W. On the Convergence Behavior of the Restarted GMRES Algorithm for Solving Nonsymmetric Linear Systems // Numerical Linear Algebra with Applications. 1994. V. 1(5). P. 427–47.
Sorensen D.C. Numerical Methods for Large Eigenvalue Problems // Acta Numerica. 2002. V. 11. P. 519–584.
Поляк Б.Т., Щербаков П.С. Робастная устойчивость и управление. М.: Наука, 2002.
Girko V.L. Circular Law // Theory Probab. Appl. 1984. №. 29. P. 694–706.
Зубов Н.Е., Рябченко В.Н. О вычислении псевдообратной матрицы. Общий случай // Вестн. МГТУ им. Н.Э. Баумана, “Естественные науки”. 2018. № 1. С. 16–25.
Зубов Н.Е., Микрин Е.А., Рябченко В.Н. О вычислении псевдообратной квадратной матрицы на основе обращения // Вест. МГТУ им. Н.Э. Баумана, “Естественные науки”. 2018. № 3. С. 24–31.
Мисриханов М.Ш., Рябченко В.Н. Матричная сигнум-функция в задачах анализа и синтеза линейных систем // АиТ. 2008. № 2. С. 26–51.
Мисриханов М.Ш., Рябченко В.Н. Размещение полюсов при управлении большой энергетической системой // АиТ. 2011. № 10. С. 129–153.
Гаджиев М.Г., Коробка В.В., Мисриханов М.Ш., Рябченко В.Н., Шаров Ю.В. Редукция уравнений электрической сети на основе матричных делителей нуля // Изв. РАН. Энергетика. 2018. № 1. С. 97–107.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Известия РАН. Теория и системы управления

Архивы выпусков Информация о журнале Отправить рукопись в журнал