Решатели математического ядра K2.SimKernel

Явный или неявный решатель ДУ моделирующей программы
Библиотека классических подпрограмм (функций), которые реализуют операцию интегрирования. (Дискретные квазианалоги интеграторов используются для решения дифференциальных уравнений).
Итерационный решатель моделирующей программы
Библиотека классических подпрограмм, которые предназначены для решения алгебраических уравнений путем итерационного подбора независимых переменных.
Оптимизирующий решатель моделирующей программы
Разновидность итерационного решателя, которая предназначена для минимизации значения функции в процессе итерационного подбора ее независимых переменных.

Математическое K2-ядро предназначено для решения систем дифференциально-алгебраических уравнений. Итерационный решатель ядра способен отыскать корни, в том числе и нелинейных алгебраических уравнений. Имеющиеся неявные дискретные квазианалоги интеграторов позволяют найти решения систем дифференциальных уравнений, в том числе и с жестким спектром, т.е. таких, которые имеют составляющие переходного процесса с существенно разными (отличающимися в тысячи раз) временными масштабами. Базовая библиотека математических блоков адаптирована к возможному динамическому изменению шага симуляции, что предоставляет пользователю возможность либо программировать шаг симуляции, либо активировать адаптивные (выбирающие автоматически величину шага) методы интегрирования. В случае решения жестких задач, когда активен итерационный решатель, эта особенность ядра позволяет существенно ускорить процесс симуляции.

Итерационный решатель математического ядра K2.SimKernel

Решить систему алгебраических уравнений [f(xn)] = 0 численно — значит найти, а точнее подогнать такие значения независимых переменных [xn], которые обнулят систему полиномов [f(xn)]. В процессе подгонки, по текущим результатам нужно принимать решение о том, какие значения искомых переменных [xn+1] подставлять в следующий раз (n – индекс текущего приближения корня). Для этого существуют разные методы.

Метод Ньютона для решения уравнения и системы уравнений:

xn+1 = xn - f(xn) / f '(xn) ,

[xn+1] = [xn] - [f '(xn)]-1 [f(xn)] .

Метод Ньютона-Рафсона предполагает возможность относительного изменения шага итерации a:

[xn+1] = [xn] - [f '(xn)]-1 [f(xn)] an .

Модифицированный метод Ньютона предполагает замораживание инверсной матрицы производных (инверсного Якобиана) на первом шаге:

[xn+1] = [xn] - [f '(x0)]-1 [f(xn)] .

Но в этом случае наблюдается лишь линейное схождение. Хотя, через несколько итераций, инверсный Якобиан можно обновлять.

Метод секущих:

xn+1 = xn - f(xn) (xn - xn-1) / (f(xn) - f(xn-1)) .

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

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

В библиотеке SimLib4Visio, для ввода свободной переменной xn в модель используется блок "неизвестная" (unknown), а блок "нулевой баланс" (=0) используется итерационным решателем K2-ядра для считывания с модели значений функции f(xn). Приведенные же выше формулы используются для расчета значения на выходе блока "неизвестная" на следующем шаге итерации. Общая форма определения системы алгебраических уравнений и требуемый для дальнейших пояснений графический материал показан на рисунке.

Итерационный процесс решения системы алгебраических уравнений методом Ньютона в K2-ядре

Алгоритм решения системы алгебраических уравнений K2-ядром, предполагает выполнение на каждом шаге симуляции следующих действий:

  1. Блоки "неизвестная" инициируются начальными значениями [xn]: либо теми которые определил пользователь, либо выбранными произвольно K2-ядром. Вызывается поток, по исполнении которого итерационный решатель считывает с блоков "нулевой баланс" матрицу столбец [f(xn)].
  2. Итерационный решатель заполняет Якобиан (матрицу производных). Для чего, уточнив по количеству обнаруженных блоков "неизвестная" порядок системы уравнений, соответствующее количество раз вызывает поток, устанавливая в очередной раз приращение лишь для одной из свободных переменных. Тем самым решатель выясняет какие свободные переменные, на какие датчики "нулевого баланса" какое влияние оказывают.
  3. Имея результаты шага 1 и шага 2, решатель находит очередное приближение итерируемых переменных [xn+1] в соответствии формулами для выбранного пользователем метода. Далее, согласно алгоритму предполагается циклический повтор описанных действий. Лишь на первом этапе блоки "неизвестная" инициируются теперь уже не начальными значениями, а найденными приближениями итерируемых переменных [xn+1]. Выход из цикла предполагается после действий первого этапа в случае, если значения матрицы столбца [f(xn)] меньше заданной пользователем ошибки итерационного подбора.

Завершая обзор итерационного решателя отметим, что  математическое ядро K2.SimKernel предназначено для моделирования динамических процессов. Это означает, что использование K2-ядра будет полноценным тогда, когда коэффициенты системы алгебраических уравнений [f(xn)] = 0 будут меняться на каждом шаге симуляции. В противном случае имеет смысл лишь один шаг и для его исполнения можно пользоваться программами-калькуляторами подобными Mathcad'у.

Решатель дифференциальных уравнений математического ядра K2.SimKernel

Дискретные квазианалоги интеграторов используются для решения дифференциальных уравнений (см. рис.). Интеграторы являются самым важным математическим блоком (1/s) в библиотеках моделирующих программ. Даже если математическое ядро не имело бы итерационного решателя, то, в целом, соответствующая моделирующая программа была бы вполне полноценным программным продуктом, который мог бы удовлетворить всем потребностям такой науки как "Теория систем автоматического регулирования"; весь математический аппарат которой опирается только на системы дифференциальных уравнений, из которых предварительно исключают все алгебраические уравнения.

Система дифференциальных уравнений в форме Коши представленная направленным графом

Широко известны и находят применение два алгоритма программ численного интегрирования: одношаговый многостадийный и многошаговый. В первом случае, на выходную величину интегратора на текущем шаге ни как не могут повлиять предыдущие значения. Второй алгоритм предполагает подобную зависимость. И тот и другой алгоритмы поддерживают несколько семейств методов интегрирования, в каждом семействе имеются методы разного порядка точности и с разными границами устойчивости. Математическое ядро K2.SimKernel поддерживает одношаговый многостадийный алгоритм (одношаговые многостадийные методы). Алгоритм предполагает разбиение установленного пользователем шага симуляции на несколько субшагов — стадий оценки предсказаний приращения интеграла ki, которые вычисляются по формуле:

(1)

ki = F (n + ci Ty[n] + j=1m∑ aij kjT .

где: F (ny[n] + ...) - функция от правой части системы дифференциальных уравнений представленных в форме Коши; T - шаг симуляции.

Итоговая же величина интеграла для n+1 точки симуляции уточняется по формуле:

(2)

y[n+1] = y[n] + i=1m∑ bi ki .

Таким образом, в действительности, математическое K2-ядро на каждом шаге симуляции вызывает поток не один раз (как видит это по работе серверов визуализации пользователь), а несколько (в зависимости от выбранного метода интегрирования).

Проанализировав формулы (1) и (2) одношагового многостадийного алгоритма, легко понять, что программа будет поддерживать альтернативные методы интегрирования, если число стадий будет другим, или будут использованы другие наборы весовых коэффициентов ci, aij, bi. Определяющие тот или иной метод интегрирования коэффициенты сводятся в таблицы Батчера, общий вид которых приведен в листинге 1.

Листинг 1 - Общий вид таблиц Батчера (John C. Butcher)

   c1  |  a11 a12 ... a1m
   c2  |  a21 a22 ... a2m
   ... |  ... ... ... ...
   cm  |  am1 am2 ... amm
 ------+-----------------
       |   b1  b2 ...  bm

 где: aij - весовые коэффициенты, учитывающие влияние предсказаний 
           приращений интеграла ki друг на друга;
      cm - масштабирующие коэффициенты для разбиения шага на субшаги;
      bi - весовые коэффициенты, используемые при вычислении итоговой
           величины интеграла для n+1 точки симуляции по вычисленным 
           предсказаниям его приращений ki;
      m - количество стадий метода интегрирования

В общем случае, уравнения (1) образуют алгебраическую систему, которую требуется балансировать, подбирая величины ki. Задача усложняется тем, что модель пользователя может так же содержать алгебраические уравнения. Другими словами, для поддержки всех возможных одношаговых методов интегрирования математическое ядро должно иметь двухкаскадный итерационный решатель. Большинство моделирующих и математических программ не имеют второго каскада решателя. Рассмотрим семейства таблиц Батчера, поддерживаемые теми или другими математическими ядрами. За каждым семейством специалисты по численным методам закрепили соответствующую абвиатуру. Перечислим те семейства, обслуживание которых в математических ядрах происходит одинаковым способом:

Отметим, что в группу DIRK методов входят так же диагональные неявные методы Рунге-Кутты с одинаковыми коэффициентами на главной диагонали (SDIRK методы). А в группу EDIRK методов входят диагональные неявные методы Рунге-Кутты с явной первой стадией и одинаковыми коэффициентами на главной диагонали (ESDIRK методы). Соответствующие списку маски таблиц Бутчера приведены в листинге 2.

Листинг 2 - Маски таблиц Бутчера

 ERK            EDIRK          DIRK           FIRK
  0 | 0  0  0    0 | 0  0  0    a | a  0  0    x | x  x  x
  x | x  0  0    x | x  a  0    x | x  b  0    x | x  x  x
  x | x  x  0    x | x  x  b    x | x  x  c    x | x  x  x
 ---+--------   ---+--------   ---+--------   ---+--------
    | x  x  x      | x  x  x      | x  x  x      | x  x  x

где: x - разные числа (не одно и тоже число)

Лишь FIRK методы вызывают затруднения в программной реализации (требуют второго каскада итерационного решателя). K2-ядро поддерживает ERK, EDIRK и DIRK методы с любым количеством стадий. Поддержка двух последних семейств делает возможным решение жестких дифференциально-алгебраических систем уравнений. При активации ERK методов K2-ядро не активирует итерационный решатель (если модель пользователя не содержит алгебраических уравнений). Для поддержки EDIRK методов K2-ядро его включает. А для DIRK-методов K2-ядро имеет специальную активируемую на один шаг симуляции версию алгоритма для итеративной установки на интеграторах начальных условий (IC).

Библиотека SimLib4Visio предоставляет пользователю возможность назначить тот или иной метод интегрирования для симуляции модели. Общий список методов можно найти в блоке SimProp (см. параметр Integration Method). Разберем вариант возможного выбора:

EDIRK43a

Здесь абвиатура EDIRK означает диагональный неявный метод Рунге-Кутты с явной первой стадией. Первая из цифр — количество стадий (4). Вторая цифра — порядок ошибки метода (3). Последняя буква введена для отличия методов, у которых первые два параметра совпадают (в подобных случаях методы могут иметь разные границы устойчивости). Для вида границы устойчивости планируется ввести третий цифровой разряд.

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

Листинг 3 - Таблицы Бутчера для поддерживаемых K2-ядром методов интегрирования

ERK11 - метод Эйлера с запазыванием
  0  |  0
-----+-----
     |  1

ERK22a
  0  |
  1  |  1  
-----+---------
     | 1/2 1/2

ERK22b - метод Рунге
  0  |
 1/2 | 1/2  
-----+---------
     |  0   1

ERK32a - метод Рунге
  0  |
 1/2 | 1/2
  1  |  0   1
-----+-------------
     | 1/4 1/2 1/4

ERK33a
  0  |
 1/2 | 1/2
  1  | -1   2
-----+-------------
     | 1/6 4/6 1/6

ERK33b - метод Хойна
  0  |
 1/3 | 1/3
 2/3 |  0  2/3
-----+-------------
     | 1/4  0  3/4

ERK44a - классический метод Рунге-Кутты
  0  |
 1/2 | 1/2
 1/2 |  0  1/2
  1  |  0   0   1
-----+-----------------
     | 1/6 2/6 2/6 1/6

ERK44b - Правило 3/8 Кутты
  0  |
 1/3 | 1/3
 2/3 |-1/3  1
  1  |  1  -1   1
-----+-----------------
     | 1/8 3/8 3/8 1/8

ERK65a - метод Кутты-Нюстрема
  0  |
 1/3 |  1/3
 2/5 | 4/25 6/25
  1  |  1/4   -3 15/4
 2/3 | (  6   90  -50  8        )/81
 4/5 | (  6   36   10  8   0    )/75
-----+-----------------------------
     | ( 23    0  125  0 -81 125)/192

EDIRK23a - метод Хаммера-Холлингсуорта
  0  |  0
 2/3 | 1/3 1/3
---------------
     | 1/4 3/4

EDIRK32a
  0  |  0
 2a  |  a   a
  1  | b1   b2  a
-------------------
     | b1   b2  a

 где: a = 1-M_SQRT2/2
      b1 = b2 = M_SQRT2/4

EDIRK43a
  0  |  0
 1/2 | 1/6 1/3
 1/2 | 1/2 -1   1
  1  |  0   0  2/3 1/3
-----+-----------------
     | 1/6 2/6 2/6 1/6

EDIRK43b - ARK3(2)4L[2]SA метод
              |
      0       |       0
              |
1767732205903 | 1767732205903  1767732205903
------------- | -------------  -------------
2027836641118 | 4055673282236  4055673282236
              |
      3       | 2746238789719  -640167445237  1767732205903
    -----     | -------------  -------------  -------------
      5       | 10658868560708 6845629431997  4055673282236
              |
              | 1471266399579  -4482444167858 11266239266428  1767732205903
      1       | -------------  -------------  --------------  --------------
              | 7840856788654  7529755066697  11593286722821  4055673282236
--------------+---------------------------------------------------------------
              | 1471266399579  -4482444167858 11266239266428  1767732205903
      b       | -------------  -------------  --------------  --------------
              | 7840856788654  7529755066697  11593286722821  4055673282236
------------------------------------------------------------------------------
              | 2756255671327  -10771552573575 9247589265047  2193209047091
      b^      | -------------  -------------  --------------  --------------
              | 12835298489170 22201958757719 10645013368117  5459859503100

EDIRK64a - ARK4(3)6L[2]SA метод
              |
      0       |       0
              |
      1       |       1              1
      -       | -------------  -------------
      2       |       4              4
              |
      83      |      8611          -1743            1
     ----     | -------------  -------------  -------------
     250      |     62500          31250            4
              |
      31      |     5012029       -654441         174375           1
     ----     | -------------  -------------  -------------  -------------
      50      |    34652500       2922500         388108           4
              |
      17      |   15267082809     -71443401      730878875       2285395            1
     ----     | -------------  -------------  --------------  --------------  --------------
      20      |  155376265600     120774400      902184768       8070912            4
              |
              |      82889           0            15625           69875           -2260             1
       1      | -------------  -------------  --------------  --------------  --------------  --------------
              |     524892           0            83664          102672            8211             4
 -------------+-----------------------------------------------------------------------------------------------
              |      82889           0            15625           69875           -2260             1
       b      | -------------  -------------  --------------  --------------  --------------  --------------
              |     524892           0            83664          102672            8211             4
 -------------------------------------------------------------------------------------------------------------
              |   4586570599         0           178811875       814220225       -3700637          61727
       b^     | -------------  -------------  --------------  --------------  --------------  --------------
              |  29645900160         0           945068544      1159782912       11593932         225920

DIRK11 - метод Эйлера с упреждением (жесткий)
  0  |  1
-----+-----
     |  1

DIRK12 - Трапециидальный метод
  0  |  1/2
-----+------
     |   1

DIRK22a L-stable
  a  |  a
  1  | 1-a  a
-----+---------
     | 1-a  a

 где: a = 1-M_SQRT2/2

DIRK33a L-stable
  a  |  a                 где: a = +.4358665215084590
  e  |  d   a                  b = +1.208496649176010
  1  |  b   c   a              c = -.6443631706844691
-----+-------------            d = +.2820667392457705
     |  b   c   a              e = +.7179332607542295

DIRK43a L-stable
 1/2 | 1/2
 2/3 | 1/6  1/2
 1/2 |-1/2  1/2  1/2
  1  | 3/2 -3/2  1/2  1/2
-----+--------------------
     | 3/2 -3/2  1/2  1/2

DIRK23 A-stable
  a  |  a
 1-a | 1-2a  a
-----+---------
     | 1/2  1/2

 где: a = (3 -/+ 3^0.5)/6

Решатели математических ядер и МегаФлопс'ики

Быстродействие динамических решателей всегда было их "больным местом". У автора нет желания развивать эту спорную тему, но умолчать о недостатке — значит скрыть информацию от пользователя. Автор оставляет за собой право написать данный параграф в свободном, не официальном стиле.

Когда неподготовленный пользователь впервые сталкивается с тем объемом вычислений, которые требуются для численного решения задач, его посещает депрессия. Рассмотрим усугубляющий оценочный пример. Пусть нам требуется смоделировать поведение нелинейной динамической системы. Какой она может иметь порядок? При решении карандашом 90% задач имеют порядок не выше третьего (если порядок выше, то проще не решать, а спроектировать систему на глазок). Такой инструмент как компьютер обязан существенно помочь специалистам. Как показывает практика, среднестатистические компьютерные модели имеют порядок от 3 до 20. Если порядок выше 20, то сам автор модели спустя месяц достаточно продолжительно разбирается в ней. Итак, у нас двадцать интеграторов. Учтем, что на один реактивный элемент обычно приходится 3..5 активных, следовательно, в дополнение к 20 дифференциальным, наша модель будет содержать около 80 алгебраических уравнений первого порядка (20*4). Теперь разберемся с тем, сколько же шагов симуляции нам нужно? Обычно получается, что 300 точек пользователи берут на сложный переходный процесс, плюс еще 600, чтобы убедится, что режим установился. Для изучения поведения моделей управляемых событиями это число увеличивается пропорционально количеству вариантов событий (состояний модели). Возьмем 10 состояний. В итоге получается 10000 шагов. Теперь накрутим процедуру оптимизации модели. Пусть мы сможем оптимизировать параметры модели за 30 симуляций. Итого: Нелинейная модель из 20 ДУ, 80 АУ, 300000 расчетных точек. Если "измерить" такую модель в блоках, то их будет около 300. На каждый блок, в среднем, в математическом ядре придется 5 операций с плавающей точкой. Т.е. поток модели будет состоять из 300*300000*5=450000000 операций с плавающей точкой. Теперь учтем, что модель придется решать, скажем, пяти-стадийным методом 450000000*5. Безусловно, на каждом субшаге нам придется заполнять якобиан размером 20+80. Это означает, что количество потоков модели увеличится в 100 раз. И, конечно, нелинейность системы алгебраических уравнений приведет к тому, что корни нам удастся найти в каждом случае, скажем, лишь за 20 итераций. Общий итог. Нам требуется выполнить 450000000*5*(20+80)*20=9.000.000.000.000 операций с плавающей точкой (+/- порядок). Такие цифры понять человеку сложно. Вот на рубеже тысячелетий фирма TI выпустила цифровой сигнальный процессор TMS320C6712 DSP. Работал он на частоте 100 МГц и имел производительность 600 MFLOPS (million floating-point operations per seconds). Получается, подобный процессор решал бы нашу задачу около 15000 секунд или 4 часа. Казалось бы, есть от чего впасть в депрессию.

Однако, в действительности, ситуация не столь удручающа. Во-первых, существуют алгоритмы построения математических ядер, которые позволяют в 3..10 раз сократить объем требуемых МегаФлопс'ов (да простит читатель автора за образный термин). Во-вторых, рассмотренная модель далека от действительной жизни. В реальном мире требуется моделировать реальные цепи преобразования энергий. Не будем говорить, что нелинейные преобразования материи неизвестны, но их явно очень мало. Положим, у вас есть два источника электрической материи (источники тока) величиной 4А  каждый. И как бы вы не выкручивались, вам ни как не составить электрическую цепь, которая бы описывалась алгебраическим уравнением второго порядка (квадратным уравнением). Другими словами из двух проводов, в каждом из которых по 4А, при заданном ограничении ну ни как не получить провода с током 4^2=16А. Такую же ситуацию мы будем иметь в любом энергетическом домене. Если соотнести сказанное с нашей моделью, то система алгебраических уравнений практически всегда будет линейной. Таким образом количество требуемых итераций мы сократили до одной, и нам теперь нужно в 100 раз ((3..10)*20) меньше МегаФлопс'ов, что делает ситуацию со временем симуляции приемлемой.

Итак, думаю, мне удалось убедить пользователя математических ядер, что фундаментальные процессы преобразования той или иной энергетической материи практически всегда линейны. Однако системы управления достаточно часто основаны на нелинейных законах. Но и это ограничение достаточно часто можно обойти. Кроме классической линеаризации, на вскидку можно вспомнить еще два метода. Во-первых, даже в случае численных решений можно пользоваться "кривыми зеркалами" — тем или иным математическим доменом (логарифмический, Лапласа, Фурье, Z-домен, ...), которые позволят заменить одну сложную операцию (нелинейную операцию) более простой (линейной). Второй способ совершенно необоснованно практически не используется. Его суть в том, что можно инкапсулировать нелинейную и линейную части модели в разных модулях и решать их разделенными математическими ядрами. К сожалению, ни одна из широко распространенных сегодня моделирующих программ не имеет разделяющегося математического ядра. Тем не менее, кроме сокращения времени симуляции, подобная техника могла бы существенно усилить идею моделирования управляемого событиями.

Резюме: Надо признать, что динамические решатели — это ненасытные пожиратели "бедных" МегаФлопс'иков :-). Однако они позволяют снизить зависимость инженера от многих математических приемов, которые используются сегодня для расчетов цепей преобразования энергий и при этом ни на ёту не приближают начинающего инженера к пониманию их физической природы.

Аспекты итерационных вызовов потока

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

Листинг 4 - Процедура вызова "потока"

  for (i=0; i < numBlock; i++) {
    MathBlock[i]->Calc();  // вызов математических функций
  }

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

Так же следует отметить особую группу блоков обладающих "эффектом памяти". В нее входят: 1/s — интеграторы, 1/z — регистры задержки, 1/e — звенья чистого запаздывания. Особенность их программной реализации в том, что соответствующие процедуры вызываются дважды. Во время первого вызова блоки устанавливают величины на своих выходах, а во время второго вызова принимают входные сигналы. K2-ядро обновляет состояние регистров и звеньев запаздывания в моменты времени кратные лишь шагу симуляции. Интеграторы обновляются либо на каждом субшаге (если активен многостадийный явный метод), либо при каждом итерационном вызове потока (если активен неявный метод). Если же блок с "эффектом памяти" (1/z+clk, 1/e+tay) имеет дополнительный вход, то увеличиваются затраты вычислительных ресурсов на идентификацию того, когда же его можно не вызывать. В результате используется более простое решение — блоки вызываются в каждом субшаге или итерации (всегда), а внутренняя процедура оборудована дополнительной проверкой счетчика стадий и флага завершения итерационного подбора. Последним заслуживающим внимания аспектом является особенность функционирования интегратора с насыщением (1/s+limit). Если этот блок находится в режиме ограничения и активен какой либо DIRK-метод, то итерационный процесс на первой стадии выполняется дважды. Сначала вычисляется баланс координат модели без ограничений, за тем, в случае выхода за пределы, алгоритм переключается, и активируется процесс итеративной подгонки выходной координаты к ограничительному пределу. На последующих стадиях мульти-стадийных методов, когда вычисляются лишь прогнозы приращений высшего порядка для величины интеграла на следующем шаге, ограничения не действуют. K2-ядро поддерживает не только параллельный, но и покаскадный вход/выход интеграторов в насыщение. При активных ERK- или EDIRK-методах интегратор с насыщением не имеет ярких особенностей. Ограничения накладываются на первой стадии.