Стабилизация плазмы на базе прогноза с устойчивым линейным приближением Текст научной статьи по специальности «Математика»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 10. 2011. Вып. 1

УДК 517.977.5+533.9.07

Е. И. Веремей, М. В. Сотникова

СТАБИЛИЗАЦИЯ ПЛАЗМЫ НА БАЗЕ ПРОГНОЗА С УСТОЙЧИВЫМ ЛИНЕЙНЫМ ПРИБЛИЖЕНИЕМ

1. Введение. Данная статья посвящена проблеме стабилизации плазмы в современных токамаках на основе управления с прогнозирующими моделями (MPC - Model Predictive Control) [1, 2]. Одним из основных преимуществ этого подхода является возможность учета ограничений на управляющие и контролируемые переменные в процессе синтеза управления на каждом такте дискретного времени [3]. Кроме того, применение MPC-подхода позволяет обеспечить высокое качество процессов управления в тех случаях, когда точная математическая модель объекта управления неизвестна. Отмеченные обстоятельства определяют обоснованность выбора управления с прогнозом для решения задач стабилизации плазмы в токамаках.

Классический вариант MPC-подхода может быть эффективно использован для формирования управления в процессе стабилизации плазмы с учетом ограничений на напряжение и силу тока в энергетической системе [4]. Тем не менее данный вариант MPC-подхода обладает также и некоторыми недостатками. Наиболее существенный из них - отсутствие гарантии устойчивости замкнутой системы. С целью устранения этого недостатка разрабатывается новый алгоритм управления, базирующийся на идеях MPC-подхода и модальной параметрической оптимизации. В рамках предлагаемого алгоритма на каждом шаге дискретного времени обеспечивается размещение собственных чисел линейного приближения замкнутой системы внутри желаемых областей комплексной плоскости, которые находятся внутри единичного круга и отражают заданные требования на степень устойчивости и колебательность замкнутой системы.

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

Эффективность предлагаемого алгоритма продемонстрирована на примере системы вертикальной стабилизации плазмы в токамаке ITER.

Веремей Евгений Игоревич — доктор физико-математических наук, профессор, заведующий кафедрой компьютерных технологий и систем факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 110. Научное направление: теория управления и ее приложения в судостроении и энергетике. E-mail: e_veremey@mail.ru.

Сотникова Маргарита Викторовна — кандидат физико-математических наук, старший преподаватель кафедры компьютерных технологий и систем факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 9. Научное направление: оптимизация динамических систем с прогнозирующими моделями. E-mail: s_margosha@mail.ru.

© Е.И. Веремей, М.В. Сотникова, 2011

2. Постановка задачи управления.

2.1. Математическая модель процесса вертикальной стабилизации плазмы в токамаке ITER. В общем виде математическая модель процесса управления плазмой может быть представлена системой дифференциальных уравнений [5]

" + “ = (1>

где Ф - вектор полоидальных потоков; R - диагональная матрица сопротивлений; I - вектор активных и пассивных токов; V - вектор напряжений, прикладываемых к обмоткам электромагнитной системы. Причем вектор Ф определяется нелинейным соотношением

Ф = Ф(1 Ip), (2)

в котором Ip - ток плазмы, а вектор выходных переменных - равенством

У = V(i,Ip)- (3)

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

Токамак ITER [6] имеет отдельный быстрый контур обратной связи для стабилизации вертикального положения плазмы. В этом контуре величина управляющего на-

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

X = Ax + bu, , . .

y = cx + du, ^ '

где x € E58 - вектор состояния; u G E1 - напряжение VS-конвертера; y G E1 - вертикальная скорость токового центра плазмы.

Необходимо отметить, что порядок рассматриваемой модели (4) очень высокий. Для упрощения решения задачи синтеза закона управления осуществим редукцию системы (4) с 58-го до 3-го порядка при помощи стандартной функции schmr пакета Matlab. В результате получим следующую передаточную функцию редуцированной SISO-модели (от входа и к выходу у):

_ 1.732 • l(T6(s - 121.1)(s + 158.2)(s + 9.641)

(s + 29.21)(s + 8.348)(s - 12.21) ' ^ >

Отметим, что полюса передаточной функции (5) определяют наиболее значимую динамику по отношению к исходной модели (4). В частности, неустойчивый полюс соответствует неустойчивости вертикального положения плазмы. Два других полюса соответствуют виртуальным токам в катушках магнитной системы. Качество процесса редукции может быть проиллюстрировано сравнением диаграмм Боде для исходной и редуцированной систем. На рис. 1 слева представлены диаграммы Боде для исходной модели (4) и редуцированной модели 3-го порядка (5), а справа - аналогичные диаграммы для исходной и редуцированной до 2-го порядка моделей. Нетрудно видеть,

что кривые для исходной и редуцированной до 3-го порядка моделей практически совпадают, в то время как они существенно разнятся для модели 2-го порядка.

10° 102 Частота, рад/с

10° 102 Частота, рад/с

Рис. 1. Диаграммы Боде для исходной (сплошная линия) и редуцированной (пунктирная) моделей

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

V™ = 0.6кУ, 7™ = 20.7Ы, (6)

где Vnax - максимальное напряжение; 1^аХ - максимальный ток в УЯ-конвертере.

Итак, линейная модель (5) совместно с ограничениями (6) будет в дальнейшем рассматриваться как базис для синтеза закона управления.

2.2. Постановка задачи оптимального управления. Будем считать, что целью управления является стабилизация вертикальной скорости токового центра плазмы. Один из возможных подходов к синтезу стабилизирующего закона управления - постановка и решение соответствующей задачи оптимизации. В рамках оптимизационного подхода задача вертикальной стабилизации плазмы может быть представлена в следующей форме. Необходимо построить такой алгоритм управления с обратной связью и = и(Ь,у), который обеспечивает минимум квадратичного функционала качества

п<Х>

.] = .] (и) = (у2 (г) + \и2 (г))йг, (7)

Jo

на движениях модели (5) при наличии ограничений (6) и гарантирует устойчивость замкнутой системы. Здесь Л - весовой множитель, определяющий соотношение между быстродействием и энергозатратами.

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

В то же время МРС-подход позволяет в явной форме принимать во внимание имеющиеся ограничения. Его базовая схема предполагает решение в режиме реального времени задачи оптимизации функционала (7) на конечном горизонте прогноза по отношению к модели (5) и ограничениям (6).

3. Управление с прогнозирующей моделью в режиме реального времени. Пусть математическая модель динамики объекта управления представлена системой нелинейных разностных уравнений в рекурсивной форме

Хі+1 = f(Хі, йі), і = к + і, і =0,1, 2,..., хи = Хк,

У = СХі . (8)

Здесь хі Є Е" - вектор состояния; йі є Ет - вектор управления; уі Є Ег - вектор выходных переменных; хи Є Е" - истинное состояние объекта управления на к-м такте дискретного времени или его оценка, полученная на базе измерений.

Рассматриваемая модель предназначена для прогноза динамики объекта управления на определенное количество шагов при заданной конечной последовательности векторов управления {й к , йи+1,..., йк+Р _1}. Подобная модель называется прогнозирующей моделью, а параметр Р - горизонтом прогноза.

В результате решения системы (8) получим конечную последовательность векторов управления {йк , йи+1,..., йк+Р_і} и соответствующую ей конечную последовательность векторов измерений {у к+1, Ук+2,..., Ук+р}. Будем говорить, что при этом сформирован прогноз движения реального объекта с горизонтом Р.

Рассмотрим функционал, характеризующий качество управления прогнозирующей моделью на горизонте прогноза:

Jk = Л (У, й), (9)

где

У = ( Ук+1 Ук+2 ... Ук+р ) Є ЕгР,

й=( йк йк+1 ... йк+р_1 )Т Є ЕтР

- вспомогательные векторы.

Известная классическая схема реализации управления с прогнозом, включенным в контур обратной связи, заключается в следующем:

1) измеряется или оценивается текущее состояние объекта Хк;

2) выбором управления на горизонте прогноза оптимизируется движение прогнозирующей модели (8) по отношению к функционалу (9);

3) найденное оптимальное управление реализуется для объекта управления только на текущем такте (от момента к до момента к + 1);

4) для следующего такта, начиная от момента времени к+1, все операции, указанные в пунктах 1)-3), повторяются заново.

По существу, в рамках приведенной схемы управление осуществляется по принципу обратной связи с дискретным поступлением информации о текущем состоянии объекта в моменты к = 0,1, 2,... .

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

Приведенная выше классическая схема допускает определенное развитие в плане использования различных вариантов реализации управления с прогнозом. Выбор

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

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

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

4. Управление с прогнозом на базе параметрической оптимизации линейного приближения. Будем считать, что математическая модель объекта управления представляется системой нелинейных разностных уравнений вида

Хк+1 = F(Х.k, йк,^к), (10)

Ук = СХк, ( )

где Ук Є Ея - вектор измеряемых переменных; Хк Є Е" - вектор состояния; йк Є Ет -вектор управления; фк Є Ег - вектор внешних возмущений.

Сформируем нелинейную прогнозирующую модель на базе уравнений (10) и представим ее следующим образом:

Хі+1 = f(Хі, йі), і = к + і, і =0,1, 2,..., Хк = Хк, ( )

Уі = СХі . (11)

Здесь Хк Є Е" - текущее состояние объекта, достигнутое на к-м такте функционирования системы. Пусть заданы последовательности векторов {г%} и {Щ}, к = 0,1,2,..., определяющие некоторое желаемое движение объекта. Осуществим линеаризацию уравнений (10) в окрестности рассматриваемого движения. В результате получим систему линейных разностных уравнений, описывающую динамику объекта в отклонениях:

Хк+1 = АХ к + Вй к + Нфк, (12)

Ук = СХ к, ( )

где Хк Є Е", йк Є Ет, Ук Є Ея, фк Є Ег - векторы состояния, управления, измерений

и внешних возмущений соответственно.

Будем формировать управление для прогнозирующей модели в виде обратной связи (регулятора) по измеряемому выходному вектору Ук:

йк = W(g, Ь)Ук, (13)

где д - оператор сдвига на такт вперед; W(д, Ь) - передаточная матрица регулятора с фиксированной структурой (заданы степени полиномов в числителях и знаменателях

всех ее компонентов); Ь € Ег - вектор настраиваемых параметров, подлежащих выбору при синтезе закона управления. Запишем уравнения прогнозирующей модели (11), замкнутой регулятором (13):

Хі+1 = f(Хі, йі), і = к + і, і =0,1, 2,..., Хк = Хк,

йі = г“ + W(q, Ь)С(Хі - гХ).

(14)

Зафиксируем некоторый вектор Ь параметров и найдем решение системы (14) для і = к, к +1,..., к + Р — 1. В результате получим последовательность векторов {Хі} (і = к +1, ...,к + Р), представляющую собой прогноз движения реального объекта с горизонтом Р.

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

который при прочих одинаковых условиях с очевидностью является функцией вектора настраиваемых параметров Ь. Здесь {х*}, г = к + 1,к + Р, и {й} г = к,к + Р — 1, -последовательности векторов состояния и управления, удовлетворяющие системе (14). Введем в рассмотрение следующую задачу параметрического синтеза:

где &и - множество настраиваемых параметров, обеспечивающих расположение корней характеристического полинома замкнутой линейной системы (12), (13) внутри заданной области Сд в единичном круге.

Следует отметить, что задача (16) является специфическим вариантом задачи нелинейного программирования со сложной целевой функцией, которая достаточно часто в практических задачах не имеет явного аналитического представления и задается алгоритмически. Кроме того, специфика задачи обусловлена наличием не менее сложных ограничений, определяемых указанными требованиями по расположению корней. Заметим, что размерность задачи оптимизации (16) зависит только от размерности вектора Ь, но не от величины горизонта прогноза Р.

Определение 1. Будем говорить, что структура регулятора (13) является полной, если степени полиномов в числителях и знаменателях компонентов передаточной матрицы W(q, Ь), а также размерность и состав компонентов вектора Ь таковы, что с помощью выбора этого вектора (т. е. назначения конкретных величин настраиваемых параметров) можно обеспечить произвольный спектр корней характеристического полинома А(г, Ь) замкнутой системы (12), (13).

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

Л = М{Хі}, {йі}) = л(W(д, Ь)) = л(Ь) > 0,

(15)

(16)

Хк+1 = Ах к + Вй к + Нфк, Ук = СХ к,

(17)

Ск+1 = Ас(Ь)£к + Вс(Ь)Ук, йк = Сс(Ь)£к + Ос (Ь)Ук,

где £к € Е" - вектор состояния регулятора. Применяя к системе (17) Z-преобразование при нулевых начальных условиях, получим

(Епг — А)х = Вй + Иср, (Е^ г — АС(Ь))£ = Вс(Ь)Сх, й = Сс(Ь)£ + Бс(Ь)Сх, У = Сх,

или

(Епг — А — ВБс(Ь)С — ВСс(Ь) \ / х\/ Н \

^ —Вс(Ь)С Е„г — АС(Ь) Д £ ) ^ 0 )*.

Отсюда следует, что характеристический полином замкнутой системы Д(г, Ь), степень которого будем обозначать н^, имеет вид

Д(г Ь) = сЫ ( Ет — А — В°с(Ь)С —ВСс(Ь)

Д(г Ь) аелу — вс(ь)с Е^г — АС(Ь)

Пусть выбором значений настраиваемых параметров, объединенных в вектор Ь, необходимо назначить заданные корни характеристического полинома системы (17), т. е. обеспечить выполнение тождества

Д(г, Ь) = Д(г),

в котором Д(г) - заданный полином степени н&. Для этого приравняем коэффициенты

при одинаковых степенях г и получим систему, состоящую из (н^ + 1) нелинейного

уравнения с г неизвестными компонентами вектора Ь:

Ь(Ь) = ъ (18)

Очевидно, что структура регулятора (13) является полной тогда и только тогда, когда система (18) совместна для любого вектора 7.

Нетрудно убедиться в том, что если в качестве компонентов вектора Ь выступают коэффициенты дробно-рациональных функций, являющихся элементами передаточной матрицы W(д, Ь) регулятора (13), то система (18) для поиска этого вектора оказывается линейной:

ЬЬ = ^, (19)

где Ь - матрица с постоянными компонентами. Как в случае (18), так и в (19) необходимыми условиями полноты структуры являются полная управляемость пары А, В и полная наблюдаемость пары А, С.

Уточним теперь постановку задачи оптимизации (16), считая, что заданная структура регулятора (13) полная. Определим допустимое множество 0,и следующим образом:

Пи = {Ь € Ег : й(Ь) € Сд, г = 1,..., ил}, (20)

здесь § - корни характеристического полинома Д(г, Ь), н^ = СедД(г, Ь).

Рассмотрим два варианта задания области Сд внутри единичного круга, представленные на рис. 2. Дадим формальное определение указанным областям:

Сд = Сд1 = {г € С1 : \г\ ^ г}, где г € (0,1) - заданное вещественное число;

Сд = Сд2 = {г € С1 : г = рехр(±гр), 0 < р < г, 0 < р < ^(р)}, где г € (0,1) -заданное вещественное число; ф(£) вещественная функция переменной £ € (0, г], принимающая значения из отрезка [0,^], причем ф(г) = 0.

1т г

\VCiZ i

/ г /Ч - ' . - "' шж

1 Щр Ке/

Рис. 2. Области Сд1 (а) и Сд2 (б) желаемого расположения корней

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

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

Рассмотрим следующее вспомогательное утверждение:

Теорема 1. Для любого вектора 7 Е ЕПа корни полинома А* (2,7), построенного по приведенным ниже формулам, находятся внутри области Сд1 или на ее границе. Обратно, если корни некоторого полинома А(г) принадлежат области Сд1 и при этом вещественные корни положительны, то можно указать такой вектор 7 е ЕПа, что справедливо тождество А(г) = А* (2,7). Здесь

А* (2,Т) = П(22 + а1 (^,г)г + а0 (7,г)),

г=1

если па - четное, й = па/2;

а

А*(2,1) = (2 - аа+1 (7,г))П(-г2 + а\(!,г)г + а0(a,r)), если па - нечетное, й = [па/2];

г=1

аг(7, Г) = -Г ( ех Р ( -^ + \/4 - 7*2

2 1 V 4 Н2 ) 1 ехР 12 V 4

о(7,г)= г2ехр'(-721) , I = !,...,й, аа+1 (7,г) = гехр (-7|0)

7«

(21)

(22)

(23)

7 = {711,712,721,722, ...,1а1,1а2,7ао}. (24)

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

трехчленов в формуле (21). Действительно, для любой пары действительных чисел 7*1, 7*2 корни г-го трехчлена в (21), с учетом (23) и (24), представляются выражением

4,2 = г ехР

4

откуда следует, что \z\2\ ^ г, т. е. корни трехчлена г\ 2 находятся внутри области СД1 или на ее границе, что и доказывает прямое утверждение.

Для доказательства обратного утверждения зададим некоторый квадратный трехчлен Д*(х) = х1 + + во. По условию будем считать, что его корни х\^ принадлежат

области СД1 и, если они вещественные, то положительны. Для того чтобы корни х1,2 принадлежали кругу \х\ ^ г, необходимо и достаточно, чтобы выполнялись соотношения [7]

(25)

Кроме того, поскольку произведение корней Х1Х2 положительно как в случае пары комплексно сопряженных корней, так и для двух вещественных положительных корней, то должно выполняться условие

во > 0.

(26)

Найдем такие 7*1 и ^и, чтобы выполнилось тождество Д*(х) = Д*(х). Приравнивая коэффициенты при одинаковых степенях, получим

в1,

откуда имеем

7*2

г2ехр(-7г21) = во,

7л = (/30/г2),

41п ГЙ) 1п

(27)

где ги = -§-о - ! + ^(^-1) -1. Убедимся

в том, что числа 7*1 и ^и, определяемые

по формулам (27), являются вещественными. Действительно, из неравенств (25) и (26) вытекает, что 0 < во/г2 ^ 1, поэтому —1п (во/г2) ^ 0 и - вещественное число.

Покажем, что подкоренное выражение для неотрицательно. Рассмотрим сначала случай, когда трехчлен Д*(х) имеет два вещественных положительных корня 2:1,2. При этом его коэффициенты должны удовлетворять условию в2 — 4во ^ 0, из которого следует, что т ^ 1 - вещественное число. В результате, с учетом (25), находим

1п (тг2/во) ^ 0.

(28)

Заметим, что из (25) вытекает также выполнение неравенства в2 — 2во ^ г2 + во/г2, при учете которого имеем

тво ^ г2, — 1п (тво/г2) ^ 0.

(29)

Из неравенств (28) и (29) следует, что подкоренное выражение для 7*2 неотрицательно и 7*2 - вещественное число.

Рассмотрим теперь случай, когда трехчлен Д*(х) имеет пару комплексно сопряженных корней Х1 ,2. При этом должно выполняться неравенство в2 — 4во < 0, из которого имеем, что т - комплексное число, представимое в виде т = в2/2во — 1 +

*причем нетрудно показать, что |гу| = 1. Тогда подкоренное выражение для 7*2 имеет вид

откуда вытекает его неотрицательность, т. е. - вещественное число.

При нечетном н^, в соответствии с (22), в состав полинома Д* дополнительно вводится линейный двучлен, для которого утверждения теоремы очевидны. ■

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

Теорема 2. Для любого вектора 7 € ЕПа корни полинома Д* (х, 7) (21),(22) принадлежат области Сд2 и обратно, если корни некоторого полинома Д(х) принадлежат области Сд2 и при этом вещественные корни положительны, то можно указать такой вектор 7 € ЕПа, что справедливо тождество Д(х) = Д* (2,7), где

ai(7, г) = —г (ехР (—12ц + ч) + ехР (—12ц — ч)) , (30)

а°(7,г)= г2ехр (—2721) ,г = 1, ...,3,, аЛ+1(^,г) = г ехр(—72о), (

где щ = ^7^ - / (7*2) (Ф2 (г ехр (-7^)) + 7а); 7 = {иъ 712, 721, 722, •••, 7<гь 7<г2, 7<го}-

При этом функция /(•) : (—те, +те) ^ (0,1) должна удовлетворять условию существования обратной функции во всей области задания, а ф(£) - вещественная функция переменной £ € (0, г], принимающая значения из отрезка [0,^], причем ф(г) = 0.

Доказательство. Как и в теореме 1, обратимся к свойствам квадратных трехчленов в (21). Вначале докажем прямое утверждение.

Для любой пары действительных чисел 7*1, 7*2 корни г-го трехчлена Д* (х) в (21) представляются выражением х\ 2 = г ехр—7^ ± V*). Здесь возможны два варианта. Если V* - вещественное число, то корни г\ 2 тоже вещественные. При этом, в силу свойств функции /, справедливо неравенство 741 — /(7*2) (ф2 (гехр(—721)) + 741) ^ 7^, из которого следует, что корни положительны и \х\ 2\ ^ г, т. е. г\ 2 € Сд2.

Если же V* - комплексное число, то г\ 2 - пара комплексно сопряженных корней, причем \х* 2 \ = р = г ехр(—7*21) ^ г. В силу свойств функции /, выполняется неравенство _______________________________ ______________

Ч>=^/(7*2) (Ф2 (г ехр(-721)) + 7^) - 7^ < ^ф2 (г ехр(-721)) = ф(р). (31)

Так как а^г\ 2 = и, с учетом (31), 0 ^ р ^ Ф(р), то корни г\ 2 принадлежат области Сд2, что и доказывает прямое утверждение.

Докажем теперь обратное утверждение. Для этого рассмотрим некоторый квадратный трехчлен Д*(х) = х2 + в1х + во, и будем считать, что его корни Х1,2 принадлежат области Сд2 , причем если они вещественные, то положительны. Заметим, что коэффициенты этого трехчлена должны удовлетворять условиям (25) и (26), поскольку \х1,2 \ ^ г и произведение корней в любом случае положительно.

Найдем такие и 7*2, чтобы выполнилось тождество Д*(х) = Д*(х). Приравнивая

коэффициенты при одинаковых степенях, получим

Покажем, что 7*2 является вещественным числом. Для рассуждения аналогичны

приведенным в теореме 1.

Очевидно, что уравнение относительно 7*2 имеет решение, если выражение, стоящее справа от знака равенства, принимает значения из интервала (0,1). Обозначим данное выражение через Н. Заметим, что знаменатель в Н обращается в нуль только в том случае, когда х1 = х2 = г, но при этом в качестве 7*2 может быть выбрано любое вещественное число. В общем случае, учитывая доказательство теоремы 1, нетрудно видеть, что Н > 0. Кроме того, справедливо неравенство

откуда вытекает, что существует вещественное число 7*2, которое является решением уравнения ](7*2) = Н.

При нечетном н^, в соответствии с (22), в состав полинома Д* дополнительно вводится линейный двучлен, для которого утверждения теоремы очевидны. ■

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

Прежде всего отметим, что при переходе от непрерывной модели к дискретной линейной системе собственные значения преобразуются по следующему правилу [8]: если в - собственное значение матрицы непрерывной системы, то х = евТ - собственное значение матрицы соответствующей ей дискретной системы, где Т - шаг дискретизации. С учетом этого обстоятельства рассмотрим примеры отображения стандартных областей на соответствующие области для дискретной системы.

Пример 1. Пусть задана область С = {в = х ± у] € С1 : х ^ —а}, представленная на рис. 3. Ясно, что точкам прямой х = —а соответствуют точки окружности \х\ = е-аТ. При этом область С отображается на круг \х\ ^ е-аТ, как показано на рис. 3, что соответствует области Сд1, определяющей заданную степень устойчивости дискретной системы.

Пример 2. Рассмотрим область, представленную на рис. 4:

где 0^/?<|иа>0 - заданные вещественные числа.

Осуществим отображение области С на плоскость г. Очевидно, что вершине угла (—а, 0) соответствует точка на плоскости г с полярными координатами г = е-аТ, ф = 0. Отобразим теперь каждый отрезок из семейства:

г (ехр (—721 + VI) + ехр (—72! — Vi)) = ві, Г2ехр(—272-і) = во,

отсюда следует, что

7*1 = Vх—°-51п(/?0/г2),

С = {в = х ± уі Є С1 : х ^ —а, 0 ^ у ^ (—х — а)tgв},

Ь1 = {в = х ± уі Є С1 : х = ^,1 < —а, 0 < у < (—7 — а^в}.

Рис. 3. Соответствие областей для непрерывной и дискретной систем (пример 1)

Рис. 4. Соответствие областей для непрерывной и дискретной систем (пример 2)

Каждой точке 5 = 7±уз отрезка Ь1 соответствует точка г = ввТ = в1Т±зуТ на плоскости г. Следовательно, точки отрезка Ь7 отображаются на дугу окружности с радиусом в'уТ, если —а — п/(Тtgв) <7 ^ —а, и на всю окружность, если 7 ^ —а — п/(Тtgв). Таким образом, наибольший радиус окружности, полностью заполняемой точками отрезка, равен т' = в-аТ-п/ь%в , соответствующий значению 70 = —а — п/(Тtgв).

Заметим, что лучи, образующие угол, отображаются в логарифмические спирали. Причем границу области на плоскости г образуют дуги этих спиралей, соответствующие изменению значения х от —а до 70. Введем обозначение р = вхТ. Определим функцию ф (р), ограничивающую значение аргумента при заданном радиусе окружности р:

ф (р) = 1 (—1пР — аТесли Р е [т' ,т],

1 п, если р Е [0, т'].

Результат отображения показан на рис. 4. Отметим, что полученная в результате отображения область определяет заданную меру устойчивости и колебательности дискретной системы.

Теперь воспользуемся теоремой 2 для построения вычислительного метода решения задачи параметрического синтеза (16) на допустимом множестве (20) при условии,

что Сд = Сд2, поскольку первый вариант с очевидностью является частным случаем

второго.

С этой целью зададим произвольный вектор 7 Е ЕПа и построим вспомогательный полином А*(г,7) по формулам (21), (22), (30). Потребуем, чтобы настраиваемые параметры регулятора (13), объединенные в вектор Ь Е Ег, обеспечивали тождество

А(г, Ь) = А* (г, 7), (32)

где А(г, Ь) - характеристический полином замкнутой системы степени ил. Приравнивая коэффициенты при одинаковых степенях г в (32), получим систему нелинейных

уравнений

ВД = хЬ) (33)

относительно неизвестных компонентов вектора Ь, которая при любых 7 Е Епл является совместной в силу полноты структуры регулятора (13). Будем считать, что в общем случае эта система имеет неединственное решение. Тогда вектор Ь может быть представлен как совокупность двух векторов Ь = {Ь, Ьс}, где Ьс Е Еп° - свободная составляющая (назначение которой произвольно), Ь - вектор, однозначно определяемый решением системы (33) при заданном векторе Ьс.

Введем обозначение для общего решения системы (33):

Ь = Ь* = {Ь *(ЬС, 7), Ьс} = Ь*(7, Ьс) = Ь*(е),

где через е = {'у, Ьс} обозначен произвольный вектор независимых параметров размерности Л, причем

Л = ёш е = ёш 7 + ёйп Ьс = + ис.

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

5ц+1 = f(х*, щ), г = к + ], ] =0,1, 2,..., х* = хк, ( )

й = г“ + W(q, Ь*(е))С(х* — г,х). (34)

При этом функционал Л* (15), вычисляемый на движениях системы (34), становится функцией вектора е:

Л* = Л* ({Х}, {й} = Л*к (W ^, Ь*(е))) = Л*(е). (35)

Тогда справедливо следующее утверждение:

Теорема 3. Если в задаче параметрического синтеза (16), где Пн - допустимое множество, определяемое соотношением (20) при условии Сд = Сд2, экстремум достигается в некоторой точке Ьк0 Е 0.н, то в пространстве ЕЛ найдется такая точка ек0, что

Ьио = Ь*(е*0), ек0 = а^ шт Л*(е). (36)

е£ЕА

И обратно: если в пространстве Ел существует точка ек0, удовлетворяющая (36), то точка Ьк0 = Ь*(ек0) является решением задачи (16). Или, иными словами, в указанном смысле задача (16) эквивалентна задаче на безусловный экстремум

Л* = Л* (е) ^ М . (37)

е£ЕА

Доказательство. Предположим, что имеет место условие

Ь*0 = а^ шт Л*(Ь), Л*0 = Лк(Ь*0). (38)

При этом замкнутая линейная система (17) будет иметь характеристический полином А(г, Ь* 0) с корнями, расположенными в области Сд2. Следовательно, по теореме 2, найдется такая точка 7 = у*0 Е Епа, что А(г, Ьк0) = А*(г,ук0), где А* - полином, формируемый по формулам (21), (22). Таким образом, в пространстве Ел существует точка е*0 = {ук0, Ьк0с}, где Ьк0с - соответствующая составляющая известного вектора Ьк0, для которой выполняются условия Ьк0 = Ь*(ек0), Л*(ек0) = Лк0.

Осталось показать, что в ЕЛ не существует точки ео1 такой, что Л*(ео1) < 3ко. Действительно, предположим обратное. Но тогда для точки Ь*(еоі) имеет место 3к(Ь*(єоі)) = 3*(ео1) < Зо, чего быть не может в силу (38). Аналогично доказывается и обратное утверждение теоремы. ■

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

Алгоритм решения задачи параметрического синтеза:

1. Задать начальную точку 7 Є ЕПа и построить полином А* (г, 7) по формулам (21), (22), (30).

2. В соответствии с тождеством А(г, Ь) = А* (г, 7) сформировать систему нелинейных уравнений

ВД= х(і), (39)

которая всегда совместна, и, если ее решение неединственное, назначить произвольный вектор свободных переменных Ьс Є Е”с.

3. При заданном векторе е = {у, Ьс} Є ЕЛ решить систему (39), получая при этом точку Ь*(е).

4. Сформировать уравнения прогнозирующей модели, замкнутой регулятором (13) с вектором параметров Ь*(е), и вычислить значение минимизируемого функционала 3*(е) (35).

5. С помощью любого допустимого численного метода решения задачи (37) на безусловный экстремум, задать новую точку е и, повторяя пункты 3, 4, минимизировать функцию 3*(е).

6. По нахождении точки еко = агё тіп 3* (е) определить вектор Ько = Ь* (еко), кото-

єЄЕА

рый и принять в качестве решения задачи (16).

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

Алгоритм управления с прогнозом на базе параметрической оптимизации:

1. Восстановить текущее состояние х к объекта с помощью асимптотического наблюдателя на основе данных измерений.

2. Решить задачу параметрического синтеза (16) в соответствии с приведенным выше алгоритмом по отношению к прогнозирующей модели (11) с начальными условиями х к = х к .

3. Регулятор (13) с найденным в результате решения задачи (16) вектором параметров Ь к о использовать на текущем такте.

4. Повторить заново все операции, указанные в пунктах 1-3, на к + 1-м такте.

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

5. Стабилизация вертикального положения плазмы на базе управления с прогнозом. Прежде всего преобразуем модель (5) к форме пространства состояний. В результате имеем

х = Ах + Ъп, (40)

у = сх + 6.П, ^ '

где х Є Е3 - вектор состояния; у Є Е1 - вертикальная скорость; п Є Е1 - напряжение в УБ-конвертере. Будем полагать, что модель (40) в точности описывает рассматриваемый процесс управления.

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

хі+1 = Айхі + Ъdпг, хк = хк, (41)

уі = CdXi. (41)

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

хк + 1 Adxk + Ъdпк, (42)

Ук = CdXk, ( )

здесь хк Є Е3, пк Є Е1, ук Є Е1. Будем формировать управление на горизонте прогноза в форме линейного пропорционального регулятора вида

пк = Кхк, (43)

где К Є Е3 - вектор параметров регулятора. Следует отметить, что в действительности управление по закону (43) осуществляется на базе оценки вектора состояния, полученной с помощью асимптотического наблюдателя. Заметим также, что регулятор (43) имеет полную структуру, так как ранг матриц управляемости и наблюдаемости для системы (42) полный.

Рассмотрим теперь уравнения прогнозирующей модели (41), замкнутой регулятором (43). В результате получим

хі+1 = А + bdK)Xi, хк = хк, (44)

Уі = CdXi. (44)

Качество процессов управления прогнозирующей моделью (44) на горизонте прогноза Р будем оценивать с помощью квадратичного функционала

р

3к = 3к(К) (Уй+3 + пік+і-1^ . (45)

3 = 1

Нетрудно видеть, что функционал (45) в действительности является функцией трех переменных - компонент вектора парметров регулятора К. Важно отметить, что минимизируемый функционал остается существенно нелинейной функцией в данном варианте МРС-подхода даже в том случае, когда прогнозирующая модель линейная. Этот факт - следствие гарантии устойчивости замкнутой линейной системы.

Рассмотрим следующую задачу параметрического синтеза:

3к = 3к(К) ^ тіп , 0>к = {К Є Е3 : 5і(К) Є Сд, і = 1, 2, 3}. (46)

КєО,к

Здесь 6і - корни характеристического полинома А(г, К) замкнутой системы (42), (43), степень которого nd = 3. В качестве области желаемого расположения корней зададим

область Сд = Сд2, где г = 0.97 и функция ф(р) имеет вид

Г 1п Г-') tg/?, если ге-/^ ^ р ^ г,

Ф(р) = 4 '

[эт, если 0 < р ^ гє-п/і‘ів,

Яег

Рис. 5. Область Сд желаемого расположения корней Пунктирной линией показан единичный круг.

где в = эт/10. Данная область (ограничена сплошной линией) приведена на рис. 5.

Сформируем систему уравнений в соответствии с тождеством А(г, К) = А*(г, 7), где 7 Є Е3 и полином А * (г, 7) находится по формулам (22), (30). В результате получим линейную систему относительно неизвестного вектора параметров К

Ьо + Ь1К = х(7), (47)

в которой вектор Ьо и квадратная матрица Ь - постоянные, полностью определяются матрицами системы (42) и не зависят от номера такта к. Кроме того, матрица Ь -невырожденная, следовательно, единственное решение системы (47) представимо в виде

К = Ьо + Ь 1х(7). (48)

Здесь Ь1 = Ь-1 и Ьо = —Ь-1Ьо. Подставляя (48) в прогнозирующую модель (41), а затем в функционал (45), получим 3к = 3к(К) = Л*(7). То есть функционал 3к сводится к функции трех независимых переменных. Тогда, в соответствии с теоремой

3, задача оптимизации (46) эквивалентна безусловной минимизации

3* = 3* (7) ^ тіп . (49)

7ЕЕ3

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

Рассмотрим теперь процессы вертикальной стабилизации плазмы.

Для начала опишем случай без ограничений на управляющие и контролируемые переменные. Напомним, что управление на горизонте прогноза формируется в виде линейного пропорционального регулятора (43). В связи с этим, если корни характеристического полинома системы (42), замкнутой LQR-регулятором, находятся внутри области Сд, то вектор параметров К будет практически совпадать с вектором коэффициентов LQR-регулятора. Корни характеристического полинома замкнутой дискретным LQR-регулятором системы следующие: г1 = 0.9591,22 = 0.8661,23 = 0.9408. Они находятся внутри области Сд. Переходные процессы в системе, замкнутой МРС-регулятором без учета ограничений, представлены на рис. 6, I.

ух,м!с

и, В

У\1 м/с

3>1,м/с

у2-10\А

II

и, В

у2- Ю4, А

III

и, В

у2-10*, А

и с

Рис. 6. Переходные процессы для МРС-регулятора без учета ограничений (I), с ограничением напряжения в УЯ-конвертере (II) и с ограничением тока и напряжения в УЯ-конвертере (III)

Рассмотрим теперь процессы стабилизации плазмы с учетом ограничений (6). Необходимо отметить, что для учета ограничения по току к системе (40) необходимо добавить дополнительное уравнение. Заметим также, что при наличии ограничений вместо безусловной минимизации (49) необходимо решать задачу нелинейного программирования с целевой функцией J(y) и ограничениями (6). На рис. 6, II показаны переходные процессы в системе, замкнутой MPC-регулятором, при учете ограничения только на напряжение в VS-конвертере. Из рисунка видно, что ограничение по напряжению выполняется, а ограничение по току нарушается. На рис. 6, III представлены аналогичные процессы для MPC-регулятора при учете ограничений одновременно на ток и на напряжение в VS-конвертере. Нетрудно видеть, что все требуемые ограничения выполняются.

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

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

Литература

1. Camacho E. F., Bordons C. Model Predictive Control. London: Springer-Verlag, 2004. 405 p.

2. Allgower F., Zheng A. (editors) Nonlinear Model Predictive Control. Basel: Birkhauser-Verlag, 2000. 472 p.

3. Maciejowski J. M. Predictive Control with Constraints. London: Prentice Hall, 2002. 331 p.

4. Sotnikova M. Plasma stabilization based on model predictive control // Intern. J. of Modern Physics A. 2009. Vol. 24, N 5. P. 999-1008.

5. Misenov B. A., Ovsyannikov D. A., Ovsyannikov A. D. et al. Analysis and synthesis of plasma stabilization systems in tokamaks // Proc. 11th IFAC Workshop. Control Applications of Optimization. New York, 2000. Vol. 1. P. 255-260.

6. Gribov Y., Albanese R., Ambrosino G. et al. ITER-FEAT scenarios and plasma position/shape control // Proc. 18th IAEA Fusion Energy Conference. Sorrento, Italy, 2000. ITERP/02.

7. Чернецкий В. И. Математическое моделирование динамических систем. Петрозаводск: Изд-во Петрозаводск. гос. ун-та, 1996. 432 c.

8. Hendricks E., Jannerup O., Sorensen P. H. Linear Systems Control: Deterministic and Stochastic Methods. Berlin: Springer-Verlag, 2008. 556 p.

Статья принята к печати 14 октября 2010 г.