Физика плазмы, 2023, T. 49, № 12, стр. 1348-1356

1. ВВЕДЕНИЕ

Для достижения желаемых конфигураций плазмы в современных токамаках необходимы системы управления положением, током и формой плазмы, также называемые системами магнитного управления плазмой. Форма плазмы не может быть измерена непосредственно и должна быть рассчитана по сигналам магнитной диагностики токамака в режиме реального времени. Эта задача называется задачей восстановления равновесия плазмы (plasma equilibrium reconstruction problem). Ранние алгоритмы восстановления равновесия имели время расчета за такт (TET) в 30 мс на токамаке DIII-D на пространственной сетке 33 × 33 [1]. Для современных алгоритмов типично TET в диапазоне от $200$ до $400$ мкс, например, 200 мкс на токамаке TCV с сеткой $28 \times 65$ [2] и 375 мкс на токамаке EAST с сеткой 129 × 129 [3]. Требования малости TET для расчета формы особенно сильны для небольших токамаков, таких как Глобус-М2 [4], с короткой длительностью разряда около $100$ мс от пробоя до вывода плазменного тока.

В этой статье описывается улучшенный алгоритм восстановления равновесия FCDI (Flux and Current Distributions Identification) [5] для работы в реальном времени на токамаке Глобус-М2 с TET в $20$ мкс, и система магнитного управления плазмой с алгоритмом FCDI в обратной связи. Система управления и алгоритм FCDI разработаны в среде Matlab/Simulink, с использованием Matlab Embedded Coder для генерации C++ кода, который компилируется в приложение реального времени для ЦМРВ Speedgoat Performance, которая может использоваться в обратной связи реальной системы управления плазмой как устройство HIL (Hardware-in-the-loop).

Система управления, представленная в этой статье, была улучшена по сравнению с [6–8], и является первой системой реального времени магнитного управления плазмой для токамака Глобус-М2 с алгоритмом восстановления в обратной связи. Эта система управления реализует согласованное управление положением и формой плазмы с автоматическим расчетом задающих воздействий на положение плазмы и обеспечивает компенсацию рассеянного поля центрального соленоида.

Объект управления описан в разд. 2. Разде-лы 3–5 описывают задачу восстановления равновесия и алгоритм FCDI. В разд. 6 показана модель объекта управления, а в разделе 7 – структура системы управления. Результаты моделирования в реальном времени показаны в разд. 8. Выводы сформулированы в разд. 9.

2. ОБЪЕКТ УПРАВЛЕНИЯ

Алгоритм восстановления равновесия плазмы и система управления, описываемые в данной статье, разработаны для использования на токамаке Глобус-М2 (ФТИ им. А.Ф. Иоффе РАН, Санкт-Петербург, Россия) [4]. Система обмоток полоидального поля токамака Глобус-М2 изображена на рис. 1 и состоит из 6 ОПП (обмоток полоидального поля): 4-секционная обмотка ОГУП (обмотка горизонтального управляющего поля), используемая для управления вертикальным положением плазмы, 2секционная ОВУП (обмотка вертикального управляющего поля), используемая для управления горизонтальным положением плазмы, ОИ (обмотка индуктора), используемая для генерации тока плазмы, 6-секционная ОК (обмотка коррекции), компенсирующая рассеянное полоидальное поле соленоида, и 2-секционные ОУ1 и ОУ3, которые предполагается использовать для управления формой плазмы. Секции обмоток соединены последовательно согласно, за исключением ОГУП, верхние и нижние секции которой соединены последовательно встречно. Магнитная диагностика токамака Глобус-М2 включает пояса Роговского, измеряющие токи в шести ОПП (${{I}_{{PF}}}$) и ток плазмы (${{I}_{P}}$), а также 21 магнитную петлю, измеряющие магнитный полоидальный поток (${{\Psi }_{{ML}}}$) вблизи камеры токамака. На данный момент в токамаке Глобус-М2 применяется аналоговая система управления положением и током плазмы и токами в обмотках, относящаяся к типу RZIP систем управления (система управления положением и током плазмы) [9].

3. ВОССТАНОВЛЕНИЕ РАВНОВЕСИЯ ПЛАЗМЫ

Токамак – аксиально-симметричная установка, поэтому равновесие плазмы в токамаке принято описывать на полоидальной плоскости $(r,z)$. Для этого вводится функция распределения полоидального потока $\psi (r,z)$ [10]. Она определяется как магнитный поток через поверхность $S$, ограниченную окружностью с координатами $(r,z)$:

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

Распределение полоидального потока $\psi $ и плотность тороидального тока $J$ в токамаке связаны дифференциальным уравнением равновесия

Плотность тока является источником в уравнении и однозначно определяет полоидальный поток, если она известна. Решение уравнения может быть найдено с использованием метода функции Грина [9]

где $K$ и $E$ – эллиптические интегралы первого и второго рода соответственно.

Таким образом, для нахождения полоидального потока следует в первую очередь найти распределение плотности тока в токамаке. Плотность тока в токамаке имеет три компоненты: плотность тока в обмотках ${{J}_{c}}$, плотность тока плазмы ${{J}_{P}}$ и плотность тока в вакуумной камере ${{J}_{{v}}}$.

Магнитная диагностика токамака Глобус-М2 включает пояса Роговского, измеряющие токи ${{I}_{c}}$ в 6 ОПП, площади поперечного сечения обмоток известны, распределение тока по сечению катушки может считаться равномерным для задачи восстановления равновесия, и следовательно, плотность тока в обмотках известна: ${{J}_{c}} = {{I}_{c}}{\text{/}}{{S}_{c}}$, $c = 1, \ldots ,6$.

Также в токамаке Глобус-М2 с помощью пояса Роговского измеряется ток плазмы ${{I}_{P}} = \iint {{J}_{P}}dS$, однако площадь плазмы не фиксирована, и распределение плотности тока плазмы неравномерно, поэтому необходимы дополнительные приближения с учетом измеренного значения тока плазмы.

Также, магнитная диагностика токамака Глобус-М2 включает 21 магнитную петлю, которые измеряют значения полоидального потока вне плазмы, накладывая дополнительные ограничения на распределение полоидального потока.

Для оценки того, насколько хорошо идентифицированные распределения $\psi $ и J = J_c + ${{J}_{{v}}}$ + J_p соответствуют измерениям, мы вводим квадратичный функционал ошибки

Здесь ${{\Psi }_{{ML{\kern 1pt} k}}}$ – это полоидальный поток, измеренный $k$-й магнитной петлей, ${{\tilde {J}}_{{v}}}$ – это оценка плотности тока в камере токамака, получение которой описывается в следующем разделе, а ${{I}_{P}}$ – это ток плазмы, измеренный поясом Роговского. Таким образом, задача восстановления равновесия заключается в нахождении распределений $\psi $ и $J = {{J}_{c}} + {{J}_{{v}}} + {{J}_{P}}$, которые минимизируют функционал (3) и удовлетворяют уравнению (1).

4. МОДЕЛИРОВАНИЕ ТОКОВ В КАМЕРЕ ТОКАМАКА

Токи в вакуумной камере (ВК) в токамаке Глобус-М2 не измеряются, однако плотность тока в ВК может быть оценена по закону Фарадея. Вакуумная камера является большой конструкцией с неравномерным распределением тока, поэтому для ее представления используются 56 элементов, каждый с собственным током ${{I}_{{v}}} = {{S}_{{v}}}{{J}_{{v}}}$, ${v}$ = 1 ,..., 56, где поперечные сечения ${{S}_{{v}}}$ элементов известны. Токи ВК возбуждаются в соответствии с дифференциальным уравнением

Здесь ${{R}_{{v}}}$ – диагональная матрица сопротивлений элементов ВК, ${{M}_{{{vv}}}}$ – матрица индуктивности для элементов ВК, ${{M}_{{{v}c}}}$ и ${{M}_{{{v}P}}}$ – матрицы взаимной индуктивности между элементами ВК и обмотками токамака и плазмой соответственно. Сопротивления и индуктивности ВК известны, как и взаимные индуктивности между ВК и обмотками, но из-за непостоянных формы и положения плазмы ее индуктивность не фиксирована и известна лишь приблизительно, данная неопределенность учитывается путем присвоения соответствующей погрешности рассчитываемым токам ${{I}_{{v}}}$.

Токи в 56 элементах ВК образуют 56-мерное пространство токов камеры. Чтобы уменьшить количество расчетов требуемых для моделирования токов камеры, в алгоритме FCDI используется разложение токов ВК на индуктивно несвязанные токовые моды. Переход к базису токовых мод должен совершаться при помощи матрицы перехода $P = \left\| {{{p}_{{ij}}}} \right\|$, ${{I}_{{v}}} = P{{\tilde {I}}_{{v}}}$, где ${{\tilde {I}}_{{v}}}$ – токи камеры в   пространстве токовых мод. Для получения матрицы перехода рассмотрим выражение тока в k-м элементе камеры через токовые моды

Поскольку искомые токовые моды индуктивно не связаны, при отсутствии внешних наводок токовые моды должны затухать экспоненциально $\tilde {I}_{{v}}^{{(j)}} \propto {{e}^{{ - t{\text{/}}{{\tau }_{j}}}}}$, где ${{\tau }_{j}}$ – характерное время затухания $j$-й токовой моды, следовательно производная тока в $k$-м элементе камеры равна

При отсутствии внешних наводок магнитный поток ${{\Psi }_{k}}$ через $k$-й элемент ВК выражается через элементы матрицы индуктивности ${{M}_{{{vv}}}} = \left\| {M_{{{vv}}}^{{(ij)}}} \right\|$

при этом его производная согласно (4) равна где $R_{{v}}^{{(k)}}$ – сопротивление $k$-го элемента камеры.

Беря производную от (7) и подставляя выражения (5), (6) и (8), получаем

Поскольку данное равенство должно выполняться при любых значениях $\tilde {I}_{{v}}^{{(j)}}$, коэффициенты перед $\tilde {I}_{{v}}^{{(j)}}$ в правой и левой части должны быть равны. Приравнивая их, получим систему алгебраических уравнений для нахождения элементов матрицы перехода $P$ и характерных времен затухания ${{\tau }_{j}}$

Несложно заметить, что данная система ставит задачу нахождения собственных векторов и значений матрицы $Q = \left\| {{{q}_{{km}}} = M_{{{vv}}}^{{(km)}}{\text{/}}R_{{v}}^{{(k)}}} \right\|$, при этом столбцы матрицы $P$ находятся как собственные векторы матрицы $Q$, а времена затухания мод ${{\tau }_{j}}$ как собственные значения матрицы $Q$.

В отличие от известных методов моделирования камеры [11, 12], использование полученного представления токов ВК через токовые моды позволяет уменьшить размерность задачи, и вместо расчета токов в 56 элементах ВК, рассчитывать только величины токовых мод с наибольшими временами затухания, показанными на рис. 2. В частности, равновесия плазмы токамака Глобус-М2, полученные с учетом 9 токовых мод с наибольшими временами затухания практически не отличаются от равновесий полученных с учетом большего числа мод, и вместо решения системы 56 дифференциальных уравнений (4) для 56 элементов ВК, достаточно решать аналогичную систему из 9-ти уравнений для 9 токовых мод ВК:

Используемые в этих уравнениях коэффициенты взаимной индуктивности между токовыми модами и обмотками токамака находятся как величины магнитного потока, создаваемые токовыми модами с единичным током. Так, согласно (5) при единичном токе в  j-й моде ток в k-м элементе ВК равен ${{p}_{{kj}}}$ и взаимная индуктивность между обмотками токамака и  j-й токовой модой равна

где $M_{{{v}c}}^{{(k)}}$ – k-я строка матрицы ${{M}_{{{v}c}}}$ из (4). Аналогично считаются индуктивности между токовыми модами камеры и плазмой. Используя выражения для магнитного потока через проводники с распределенным током, получим матрицу взаимной индуктивности между токовыми модами. Взаимная индуктивность между k-й и  j-й токовой модой равна подставляя (9), получим но поскольку матрица индуктивности симметрична $M_{{{vv}}}^{{(mn)}} = M_{{{vv}}}^{{(nm)}}$, можно также записать

В случае $k \ne j$, в общем случае ${{\tau }_{j}} \ne {{\tau }_{k}}$, и из (11) и (12) следует $\tilde {M}_{{{vv}}}^{{(kj)}} = 0$, т.е. полученные моды действительно индуктивно не связаны. В случае $j = k$ мы получаем выражение для собственной индуктивности k-й моды

При этом эффективное сопротивление $k$-й моды равно

5. АЛГОРИТМ FCDI-FF

Плотность тока в обмотках ${{J}_{c}}$ известна из измерений, а токи в ВК могут быть оценены согласно (10). Поэтому основной задачей восстановления равновесия плазмы является определение плотности тока плазмы ${{J}_{P}}$. Чтобы удовлетворять требованиям к скорости работы алгоритма для токамака Глобус-М2, был разработан режим FF (Fixed Filaments) алгоритма FCDI [5], реализующий метод токовых колец (филаментов) [13, 14] для восстановления равновесия плазмы. В этом методе плазма аппроксимируется рядом токовых колец с заданными пользователем положениями, распределение плазменного тока и порождаемого им полоидального потока при этом выглядят следующим образом:

Токи в филаментах определяются минимизацией функционала ошибки (3).

Так как восстановление равновесия это обратная задача, она поставлена некорректно по Адамару, и минимизация функционала (3) требует регуляризации для фильтрации физически нереалистичных решений с большой нормой. Для этой цели алгоритм FCDI использует метод SVDT (SVD Truncation) [15]. В этом методе квадратичный функционал (3) записывается в матричной форме ${{\left\| {Ax - b} \right\|}^{2}}$, где $x$ – это вектор свободных параметров (токи ВК и филаментов). Решения с большими токами фильтруются при отбрасывании вкладов малых сингулярных значений матрицы $A$ в решение.

Результаты работы алгоритма включают координаты 24 точек на границе плазмы (с шагом 15° вдоль окружности), расстояния между этими точками и лимитером (первой стенкой ВК) токамака, вытянутость и треугольность плазмы, а также положения двух ударных точек ${{g}_{{1 - 2}}}$ (точек пересечения лимитера и ограничевающей плазму линии уровня потока, называемой сепаратрисой, см. рис. 1).

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

Алгоритм FCDI-FF был внедрен на ЦВРМ Speedgoat Performance с процессором Intel i7-7700K, работающим на частоте 4.2 ГГц. Результаты моделирования в реальном времени показали TET в 20 мкс при использовании сетки из 1000 точек и 14-ю филаментами, что достаточно быстро для управления формой плазмы в токамаке Глобус-М2.

6. ДИСКРЕТНАЯ LPV-МОДЕЛЬ ПЛАЗМЫ

Восстановленные распределения плотности тока плазмы ${{J}_{P}}$ использовались для расчета массива линейных моделей ${{\{ A,B,C\} }_{k}}$, описывающих динамику плазмы в разряде. Здесь индекс $k$ обозначает момент времени ${{t}_{k}} = {{T}_{s}}k$, для которого получается модель [16]. Массив линейных моделей образует LPV-модель (Linear Parameter Varying)

Вектор входов $u$ состоит из отклонений напряжений на обмотках $u = \delta {{U}_{{PF}}}$. Вектор выходов $y$ включает в себя вертикальное и горизонтальное смещения плазмы, тока плазмы, отклонения токов в шести обмотках и значений полоидального магнитного потока на 21-й магнитной петле, $y = \,{{\left[ {\delta Z\,\delta R\,\delta {{I}_{P}}\,\delta I_{{PF}}^{T}\,\delta \Psi _{{ML}}^{T}} \right]}^{T}}$, σI_PF = [δI_HFC δI_VFC δI_CS $\delta {{I}_{{CC}}}\,\delta {{I}_{{P{{F}_{1}}}}}\,\delta {{I}_{{P{{F}_{3}}}}}{{]}^{T}}$. Состояние системы описывается вектором $x$. Шаг дискретизации модели ${{T}_{s}} = $ = 100 мкс. Всего LPV-модель имеет 30 выходов, 6 входов и 24 состояния.

7. СТРУКТУРА СИСТЕМЫ УПРАВЛЕНИЯ

На рис. 3 показана блок-схема иерархической каскадной системы магнитного управления плазмой в токамаке Глобус-М2. Внешний каскад управления содержит многомерный матричный ПИД-регулятор ${{C}_{{g\& {{I}_{p}}}}}$ для управления током и формой плазмы. Выходами этого каскада являются задающие воздействия на вертикальное $\delta {{Z}_{{{\text{ref}}}}}$ и горизонтальное $\delta {{R}_{{{\text{ref}}}}}$ положения плазмы, что обеспечивает согласованное управление положением и формой плазмы. Также на выходе этого каскада находится вектор задающих воздействий на отклонения токов в обмотках $\delta {{\tilde {I}}_{{P{{F}_{{{\text{ref}}}}}}}}$ = = ${{\left[ {\delta {{I}_{{C{{S}_{{{\text{ref}}}}}}}}\;\delta {{I}_{{P{{F}_{{{\text{1ref}}}}}}}}\;\delta {{I}_{{P{{F}_{{{\text{3ref}}}}}}}}} \right]}^{T}}.$ Вектор ${{\left[ {{{I}_{{{{P}_{{sc}}}}}}\;I_{{P{{F}_{{sc}}}}}^{T}\;\Psi _{{M{{L}_{{sc}}}}}^{T}} \right]}^{T}}$ содержит сценарные значения тока плазмы, токов в 6 обмотках и потоков на 21 магнитной петле. ${{g}_{{sc}}}$ – сценарные значения двух точек ударных точек и 3 зазоров между плазмой и лимитером (см. рис. 1), выбранные для управления из 24 зазоров, рассчитываемых алгоритмом FCDI-FF. Внутренний каскад управления содержит два одномерных ПИД-регулятора ${{C}_{Z}}$ и ${{C}_{R}}$ для управления положением плазмы и один многомерный матричный ПИД-регулятор ${{C}_{{PF}}}$ для управления токами в обмотках.

На рис. 4 показана блок-схема внутреннего каскада управления с источниками питания обмоток. Обмотка коррекции (ОК) должна компенсировать рассеянное полоидальное поле, создаваемое обмоткой индуктора (ОИ). Синтезированная система управления обеспечивает пропорциональность токов в ОИ и ОК, т.е. ${{I}_{{CC}}} \approx k{{I}_{{CS}}}$, где коэффициент $k$ зависит от сценария. Исполнительными устройствами являются источники питания обмоток. Блоки ${{W}_{{HFC}}}$ и ${{W}_{{VFC}}}$ на рис. 4 обозначают линейные модели автоколебательных инверторов тока, а блоки ${{A}_{{CS}}}$, ${{A}_{{CC}}}$, ${{A}_{{P{{F}_{1}}}}}$ и ${{A}_{{P{{F}_{{\text{3}}}}}}}$ обозначают линейные модели тиристорных выпрямителей.

Все ПИД-регуляторы в системе управления были синтезированы методом линейных матричных неравенств с использованием подхода [17]. Эффективность таких регуляторов в непрерывном времени уже была продемонстрирована для иерархической системы управления положением, током и формой плазмы в токамаке Глобус-М2 [8] и для каскадной системы управления положением плазмы в токамаке IGNITOR [18]. Этот метод был модифицирован для выполнения синтеза в дискретном времени, который ранее использовался для системы управления положением плазмы в токамаке T-15МД [19]. Метод используется для синтеза многомерных матричных ПИД-регуляторов с $m$ входами и $q$ выходами, имеющих дискретную передаточную функцию

где коэффициенты $\left\{ {{{K}_{P}},{{K}_{I}},{{K}_{D}}} \right\} \in {{\mathbb{R}}^{{m \times q}}}$ являются недиагональными (произвольными) матрицами, а $z$ – переменная $Z$-преобразования.

Этот метод не требует предварительной развязки каналов управления, обеспечивает настройку системы управления в определенной частотной области и позволяет задавать ограничения на ${{\mathcal{H}}_{\infty }}$-нормы различных Передаточных функций замкнутой системы. Сначала синтезируется регулятор ${{C}_{{PF}}}(z) \in {{\mathbb{C}}^{{6 \times 6}}}$, затем регуляторы ${{C}_{Z}}(z) \in {{\mathbb{C}}^{{1 \times 1}}}$ и ${{C}_{R}}(z) \in {{\mathbb{C}}^{{1 \times 1}}}$, и наконец регулятор ${{C}_{{g\& {{I}_{P}}}}}(z) \in {{\mathbb{C}}^{{6 \times 5}}}$. Каждый регулятор синтезируется на LPV-модели плазмы (15). В результате получается один набор регуляторов, который обеспечивает достаточный запас устойчивости и требуемое качество управления для каждой модели из массива ${{\{ A,B,C\} }_{k}}$.

8. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ В РЕАЛЬНОМ ВРЕМЕНИ

Моделирование проводилось на стенде реального времени [20] состоящим из двух ЦМРВ Speedgoat Performance, соединенных в обратную связь “Модель объекта $ \rightleftarrows $ Регулятор” через АЦП и ЦАП. На рис. 5 показана блок-схема при моделировании системы уравления в реальном времени.

Согласно данной концепции, внутри блока ЦМРВ “Регулятор” находится внешний каскад управления током и формой плазмы с алгоритмом восстановления. Это позволяет использовать новую цифровую систему управления током и формой плазмы в экспериментах, где ЦМРВ “Модель объекта” заменяется реальным токамаком вместе с существующей аналоговой системой управления положением и током плазмы.

Другая концепция предлагает размещение обоих каскадов управления в блоке ЦМРВ “Регулятор”, при этом существующая аналоговая система управления положением и током плазмы заменяется новой цифровой системой управления.

На рис. 6 показаны результаты моделирования в реальном времени, где система управления смещает ударные точки ${{g}_{{1 - 2}}}$ и зазоры ${{g}_{{3 - 5}}}$ из одного положения в другое во время диверторной фазы плазменного разряда #37326. Система управления обеспечивает слежение за задающими воздействиями ${{g}_{{1 - 6\;{\text{ref}}}}}$ и удовлетворительное качество управления.

9. ЗАКЛЮЧЕНИЕ

Моделирование в реальном времени подтвердило работоспособность разработанной системы магнитного управления плазмой, поэтому она готова к использованию в реальном эксперименте. Максимальное TET составляет $78$ мкс, а среднее TET $56$ мкс. Период дискретизации ${{T}_{s}} = 100$ мкс. Примерно $20$ мкс тратится на работу алгоритма восстановления равновесия, около $5$ мкс – на регуляторы, а оставшиеся $42$ мкс уходят на преобразование сигналов в АЦП и ЦАП. Система управления удовлетворяет условию работы в реальном времени ${\text{TET}} \leqslant {{T}_{s}}$.

Для использования такой системы управления в реальном эксперименте необходимо выполнить следующую процедуру:

– рассчитать сценарий и провести плазменный разряд,

– восстановить равновесие плазмы и рассчитать дискретную LPV-модель плазмы,

– синтезировать систему управления на этой модели и выполнить моделирование в реальном времени, чтобы убедиться в ее работоспособности и качестве управления,

– применить разработанную систему магнитного управления плазмой с алгоритмом восстановления равновесия в эксперименте.

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

Работа выполнена при финансовой поддержке РНФ (проект № 21-79-20180). Экспериментальная часть была выполнена на УНУ “Сферический токамак Глобус-М”, входящей в состав ФЦКП “Материаловедение и диагностика в передовых технологиях” ФТИ им. А.Ф. Иоффе.