ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2013 Управление, вычислительная техника и информатика № 3(24)
УДК 519.2
М.Ю. Приступа, В.И. Смагин
ПРОГНОЗИРУЮЩЕЕ УПРАВЛЕНИЕ ВЫХОДОМ НЕСТАЦИОНАРНОЙ ДИСКРЕТНОЙ СИСТЕМЫ ПРИ ОГРАНИЧЕНИЯХ НА УПРАВЛЕНИЕ1
Рассматривается задача синтеза прогнозирующего управления, построенного на основе слежения за выходом системы для нестационарного линейного объекта с учетом ограничений на управление. Прогнозирование осуществляется на основе оценок состояний нестационарного объекта, построенных с использованием экстраполятора Калмана и оценок неизвестного входа.
Ключевые слова: дискретные нестационарные системы, прогнозирующее управление, управление выходом.
При синтезе управлений широко используется метод управления динамическими объектами с применением прогнозирующих моделей - Model Predictive Control (MPC) [1, 2]. Область применения метода MPC охватывает задачи управления производственными системами, управление запасами и финансовую математику [3-8].
В данной работе рассмотрено прогнозирующее управление нестационарным объектом, модель поведения которого описывается линейными разностными уравнениями. На объект наложены ограничения, которые представлены в виде неравенств. Синтез управления с прогнозирующей моделью осуществляется на основе управления выходом объекта. Целевая функция, представленная через выпуклую квадратичную функцию, предполагает отслеживание заданного сигнала и пересчитывается в каждый текущий момент времени. Настоящая статья является продолжением работ авторов [7-10].
1. Постановка задачи
Пусть модель нестационарного объекта, канала наблюдений и управляемого выхода описываются следующими соотношениями:
xt+i = Atxt + Btut + wt, xt\t=o = xo ; ()
Vt = Htxt + vt ; ()
yt = Gtxt, (3)
где xt e Rn - состояние объекта, ut e Rm - управляющее воздействие (известный
вход), e Rl - наблюдения, выход системы контроля, yt e Rp - управляемый
выход, At, Bt, Ht, Gt - матрицы соответствующих размерностей. Уравнение (2)
является моделью системы контроля за состоянием объекта.
Далее будем полагать, что случайные возмущения wt и шумы измерения vt подчиняются гауссовскому распределению с нулевым средним и с соответствую-
1 Работа выполнена в рамках государственного задания Минобрнауки РФ на проведение научных исследований в Томском государственном университете на 2012-2014 годы, задание 8.4055.2011.
щими ковариациями W, V, то есть
Ы{у,} = 0, Ы{у,} = 0 , Ы{уут } = WtStJt, Ы{у,у1 } = V(5а , Ы{уТ } = 0. (4)
В (1) вектор начальных условий Хо является случайным, некоррелированным с величинами у, и V, и определяется следующими характеристиками:
Ы{Хо} = Хо, Ы{(Хо -Хо)(Хо -Хо)т} = РХг . (5)
Ограничения на управление представляются в виде следующих неравенств:
а1 (,) < Sut < а2 (,), (6)
где S - структурная матрица полного ранга, состоящая из нулей и единиц, определяющая компоненты вектора и,, на которые накладываются ограничения; а1 (,),
а2 (,) - заданные вектор-функции соответствующих размерностей.
Кроме того, предполагается, что пара матриц А,, В, управляема, а пара матриц А,, П г полностью наблюдаема.
Модель (1)-(3) реализует прогноз поведения объекта на некоторый период, который называется горизонтом прогнозирования и обозначается Ы, используя информацию об управлении и, и векторе наблюдений у, до текущего момента времени .
Задача состоит в том, чтобы по наблюдениям у, определить стратегию управления, при которой вектор выхода системы у, будет близок к заданному вектору у, с учетом ограничений (6).
2. Прогнозирующая модель
Поскольку случайные возмущения у и шумы измерения V, имеют гауссовское распределение, то можно выполнить оптимальное прогнозирование поведения
объекта и вектора выхода, используя экстраполятор Калмана [11]:
хі+1|, = Л%-1 + В,и, + К, (у, - ПЛ|,-1), Хо|-1 = Хо; (7)
у,+1|, = °л+1|,; (8)
К = А,РП,т (П PH, т+ V,)-1; (9)
Р+1 = Wt + АРА1-АРП,т (П,РП,Т+ V, )-1 ПРА , Ро = Рхо, (1о)
где и Уг+\\г — оценки состояния и вектора выхода в момент времени г+1,
уравнение для Р г (10) известно как разностное уравнение Риккати с дискретным временем, Рх — начальное значение дисперсионной матрицы.
Для реализации модели прогнозирующего управления необходимо иметь возможность вычислять оценки вектора состояния на моменты времени г+1, г+2, ..., г+Ы, основываясь на информации, имеющейся в момент времени г. Из уравнений (7)—(10) можно получить хг+1|г, а также оптимальные оценки для моментов
г+2,..., г+Ы:
х+М\г = А+1 +\ + вмЩ+ц , -х1|0 = х0 , I = 1Ы -1; (11)
Угщ=6м хг, * = 1 Ы . (12)
Здесь обозначение им^ используется для того, чтобы отличать действующее управление в момент г+г им от тех, которые используются с целью прогнозирования — Щ+\г.
Управляющие воздействия, используемые с целью прогнозирования, ищутся на горизонте управления М, а прогнозирование состояния объекта проводится на большем промежутке N (М < N ). Для расчета хг+м+2\г,., х+щ в (11) в качестве иг+м+1\г,..., и+Ы^ используется и1+М\г и держится на достигнутом постоянном уровне.
Уравнение (11) может быть записано через начальное состояние хг+1 \г и будущие управляющие воздействия им \ г следующим образом:
хг+2\г = Аг+1 хг+1|г + Вг+1иг+1|г ,
хг+3|г = Аг+2 хг+ 2|г + Вг+2 иг+2|г = Аг+2(Аг+1 хг+1|г + Вг+1иг+1|г) + Вг+2иг+2|г =
= Аг+ 2Аг+1 хг+1|г + А+2Вг+1иг+1|г + Вг+2иг+2|г,
i-1
=|П A
A't+i\t | 11 t+i-k I t+1\t
кк=1 ) к=1V l=1
+ Х | П А+i-i I-kuM-k\t , i = 2, Ж . (13)
Аналогично может быть преобразовано уравнение (12):
yt+1 \t = Gt+1 xt+1 \t,
yt+2\t = Gt +2 Xt +2\t = Gt+2 (At+1 xt+1\t + Bt+1ut+1\t ) = Gt+ 2At+1 xt+1\t + Gt+ 2Bt+1ut+1\t , •yt+3\t = Gt +3Xt+3\t = Gt+3 (At +2 At+1 ■Xt+1\t + At+ 2Bt+1ut+1\t + Bt +2ut+ 2\t ) =
= Gt+3 At+ 2 At+1Xt+1\t + Gt+3 At+ 2 Bt+1ut+1\t + Gt+3 Bt+ 2ut+2\t = ,
yt+i\t Gt+i | 1 1 At+i-k | Xt+1\t Gt +i | || At+i-l | Bt
к-1
, Bt+i-kut+i-k|t , г' = 2, \ . (14)
vk=1 ) k=1 V l=1
Уравнения для прогнозируемых векторов состояния и выхода могут быть представлены в векторно-матричном виде. Для этого введем следующие обозначения:
XXt =
Xt+1!t " " yt+1\t " ut+1\t
, Yt = , Ut =
Xt+N\t lyt+N\t _ ut +M\t _
En
A,
t+1 At+2 At+1
N-1
П A
L k=1
t+N-k
л, =
Gt+1
Gt+2 At+1 Gt+3 At + 2 At+1
N-1
G+N П At
t+N-k
k=1
(15)
В случае, когда горизонт управления не равен горизонту прогнозирования (М < N ), матрицы Рг и Фг вводятся следующим образом:
Рг =
0
Б,
'г+1
'г+1
Б
'м-1 \ (м _2
П А +м_к+1 I Бг+1 1П А-к=1 ) V к=1
м \ Гм-1
ПА+м_к+2 I Бг+1 |П Аг к=1 ) V к=1
П А-к=1
^+N_k |Бг+1
'г+2
+м_к+1 | Бг+2 '
г+м_к+1 | Бг+2
П А
к=1
+Ы_к | Бг+2
Б+
А+м+1Бг+м
(Ы_м _1 р
X ПАг+ы_к р к=1
Б
'г+м
Фг =
'г+2 Бг+1 'г+3 Аг+2 Бг+1
'м _1
к=1
0
0
'г+3Бг+2
П А+м_к | Бг+1 'г+м+1
'м _2
П А-
к=1
^ м \ [м _1
+м _к | Бг+2 '
'г+м+21 ПАг+м_к |Бг+1 'г+м+21 П Аг+м_к |Бг+ V к=1 ) V к=1
ГЫ_3
&г+ы |П Аг+Ы_к |Бг+1 &г+ы |П Аг+Ы_к | Бг+2 Gг+N X ПАг+ы_к
ч к=1 ) V к=1 ) V р=1 к=1
'г+м+1Бг+м
''г+м+2 Аг+м+1Бг+м
(Ы_м _1
Б
'г+м
. (16)
В случае равенства горизонта управления горизонту прогнозирования (м = N) в качестве Рг и Фг берутся матрицы следующего вида:
р, =
0
Бг+1
А + 2 Бг+1
N _2
П А
IV к=1
г+N _ к | Бг+1
N _3
П А
к=1
г+N_к | Бг+2
Б
г+N _1
0
'г+2 Бг+1 'г+3 А + 2 Бг+1
к=1
0
0
'г+3 Бг+2
П Аг+N _к | Бг+1 'г+N
П А
к=1
-г+N_к | Бг + 2
'г+NБг+N _1
Таким образом, прогнозирующая модель опишется следующей системой:
*і = * Л+ш+ Ри;
^г= ЛА+ш + Фи ■
(18)
(19)
Ограничение (6) может быть также преобразовано в векторно-матричном виде:
$1 (г +1) 5иг+1|г ^2 (г +1)
:<:<:. (20) г1(г + N)_ ^мг+м-|г _ а2(г + N) _
Введем обозначения:
а1 (і +1) а2 (і +1)
« (і) = , «2 (і) =
а1 (і + N) а2 (і + N)
5 = ^(5, ..., 5).
м
Тогда ограничения (20) для прогнозирующей модели (18)-(19) запишется следующим образом:
аг(г)<5и( <а2(г). (21)
3. Синтез прогнозирующего управления
В качестве целевой функции выберем критерий следующего вида:
N л м
' к=1
• к=1
иі+к| і иі+к-1|і||д
(-- )
где уі+к - вектор значений отслеживаемого сигнала в момент времени і + к, Сі и Бі - симметричные положительно определенные матрицы.
Отметим, что критерий (22) является классическим в задачах синтеза прогнозирующего управления [1, 2]. Квадратичная форма
IIу - у ІІС, =(у- у)Т с і(у і- у і)
обеспечивает механизм, допускающий взвешивания выходов.
Аналогично расписывается квадратичная форма ||иг+к,г _ иг+к^ :
||иг+к|г _ иг+к_1|г||0( = (иг+к|г _ иг+к_1|г) Аг (иг+к|г _ иг+к_1|г),
что предусматривает штрафы для управлений.
В случае, когда желаемая отслеживаемая траектория уг+к не известна для к > 0, целесообразно считать уг+к = у, то есть предполагать, что заданный уровень держится на протяжении всего времени управления постоянным.
Введем вектор Уг:
У+1 _ Уі+ы _
Тогда первое слагаемое в (22) преобразуется следующим образом:
1 ^ II у — ||2 1 ||у у\|2
7х\\Уг+к\г _ —г+к\\с = 2Гг _¥4с =
2 к=1 ' 2 '
= 2иТфТСгФи + и [фТС,Л,Хг+1|г _ФТС,^] + «1 .
(24)
Здесь аг - постоянная составляющая, которая не зависит от П( и хг+111, а Сг имеет
вид
[С, 0 0 ]
Сг = 0 Сг 0 . (25)
_ 0 0 С, _
Вторая часть критерия (22) может быть преобразована следующим образом:
IV.1
XII
к=1
иг+к |г иг+к _1|г|А
= ||иг+1|г иг| г|А +^ °|, иг
лг+ 2|г мг+1|г
+... + и
г +м|г иг+м _1|г
= им\АЩ+1|г _ 2и^^г|гDгuг + и1 °ги1 +
+uг+2|гDгuг+2|г _ ut+2|гDгuг+1|г + иг+1|гАгиг+1|г _ иг+1|гАгиг + 2|г + +uг+3|гDгuг+3|г _ иг+3|г°гиг+2|г + иг+2|гDгuг + 2|г _ иг+2|гDгuг+3|г +
т ^ ^ ^
+иг +м\гАгиг+м|г _ иг+м\гАгиг+м_1|г + иг+м_г|гDгuг+м_1|г _ иг+м_г|гDгuг+м|г.
Путем группировки слагаемых полученное выражение может быть записано в виде
IV.1
XI
к=1
г+к\г иг+к_1\г
А = и, 0,4 _ 2иг+1|гАгиг +
+2иг+1|гАгиг+1|г иг+1|г°гиг+2|г '
иг+2|гDгuг+1|г + 2uг+2|гDгuг+2|г иг+2|г°гиг+3|г +
иг +м _г\гDгuг+м _2|г + 2иг +м _г|г0гuг+м _1|г иг +м _1| г°гиг +м|г
_и,
-г\tDtu^
г +м\г + 2иг+м\г0гиг +м\г •
В векторно-матричном виде второе слагаемое квадратичной функции может быть записано так:
0,^11иг+к|г иг+к_1|г||0 ~ иг Агиг иг+1|г0гиг + а ;
2 к=1 А 2
где а ( - постоянная, не зависящая от и (+щ (к = 1, М ); А - блочная матрица вида
0 0
А =
2 А -А о
-А 2Аі -А
-А 2А -А 0 -А 2а
(27)
Таким образом, с учетом сделанных преобразований целевая функция запишется следующим образом:
(28)
где аі есть линейная комбинация а! и а?,
Р =ФТ С, ф+А, /і =Г
Аи
0
Гі=[фТС,Лі -ФТСі ] ■ (29)
Аналитическое решение данной задачи без учета ограничений можно полу чить из условия
& ( +цг,и г)
и
= 0
с использованием формул векторно-матричного дифференцирования [12]
т т ёХгАХВ . т ё\хАХ В
у Ау = ііАуу
ёХ
= АтВт
ёХ
= ВА.
(30)
(31)
В результате получится:
д1 (хі+і\і,иі) = д ди = ди
1 иТри + иТ/ + Сі
1 д( ігри и) +_д (иТ/)
ди
ди
= 1 [ Ртиі + ри ]+/ = 0.
(32)
В силу симметричности матрицы уравнение (32) перепишется в виде
+ А = 0.
Аналитическое решение этого уравнения имеет вид
и* = -р-1/ = -(Ф,ТСіФі + А)-1(Ф,ТСіЛХі+Ці -ФТСД)-
К 0 у
(33)
следовательно, оптимальное прогнозирующее управление записывается в следующем виде:
и+1іі= (Еп 0 ... 0)и(.
(34)
Учет ограничения (20) может быть выполнен численно, например, используя для оптимизации (28) процедуру quadprog системы МаНаЪ.
х,+1 = 0,85 + 0,05зт(0,8/) 0 _ 0,8 0,95_ х + " 0,15" 0,25_ и , х0 = "0" _0_
У, = "1 0" _0 1_ х, +V
4. Моделирование
В данном разделе приведен пример синтеза управления нестационарным объектом второго порядка с прогнозирующей моделью на основе слежения за выходом:
(35)
(36)
* =[1 0] х,, (37)
где х, є Я2 - вектор состояния, и, є Я1 - управление, у, є Я2 - вектор наблюдений, у, є Я1 - вектор выхода, шумы измерения V, являются гауссовскими величинами с нулевым средним и ковариационной матрицей
"0,0005 0 "
0 0,0005
V =
В данном примере ограничения на управление представлены следующими неравенствами:
-0,5 < щ < 0,5. (38)
Моделирование проведено на временном промежутке в 100 единиц в предположении равенства горизонтов управления и прогнозирования (М = N = 10).
Цель моделирования - определить такую последовательность управления, при которой будет отслеживаться следующая траектория:
Уи =
(39)
[-0,3 при t = 31,70
В качестве весовых коэффициентов выхода объекта и управления (в данном случае используются именно коэффициенты, а не матрицы, поскольку размерно-
сти векторов выхода у, и управления
равны единице, то есть
Шш( у.) = Шш(щ) = 1) использованы С = Д = 1.
На рис. 1 представлена динамика поведения первой (отслеживаемой) компоненты вектора состояния.
Рис. 1. Динамика первой компоненты состояния:
1 - заданная траектория (желаемая), 2 - первая компонента состояния
и
На рис. 2 приведен график управляющих воздействий на объект.
Рис. 2. Динамика управления объектом:
1 - ограничения на управление, 2 - управление
Из рис. 1 видно, что желаемая траектория отслеживается, при этом ограничения на управление соблюдаются (рис. 2).
Заключение
Разработаны алгоритмы синтеза прогнозирующего управления выходом дискретных нестационарных объектов со случайными возмущениями в условиях косвенных наблюдений за состоянием при ограничениях на управляющие воздействия.
ЛИТЕРАТУРА
1. Camacho E.F., Bordons C. Model predictive control. London: Springer-Verlag, 2004. 405 p.
2. Maciejowski J.M. Predictive control with constraints. Prentice Hall, 2002. 331 p.
3. Перепелкин Е.А. Прогнозирующее управление экономической системой производства, хранения и поставок товара потребителям // Экономика и математические методы. 2004. Т. 40. № 1. С. 125-128
4. Домбровский В.В., Домбровский Д.В., Ляшенко Е.А. Управление с прогнозирующей моделью системами со случайными зависимыми параметрами при ограничениях и применение к оптимизации инвестиционного портфеля // Автоматика и телемеханика. 2006. № 12. C. 71-85.
5. Смагин В.И., Смагин С.В. Управление запасами по двум критериям с учетом ограничений // Вестник Томского государственного университета. 2006. № 290. С. 244-246.
6. Смагин В.И., Смагин С.В. Адаптивное управление запасами с учетом ограничений и транспортных запаздываний // Вестник Томского государственного университета. Управление, вычислительная техника и информатика. № 3(4). 2008. С. 19-26.
7. Киселева М.Ю., Смагин В.И. Управление производством, хранением и поставками товаров на основе прогнозирующей модели выхода системы // Вестник Томского государственного университета. Управление, вычислительная техника и информатика. 2009. № 2(7). С. 24-31.
8. Приступа М.Ю., Смагин В.И. Прогнозирующее управление дискретными системами с неизвестным входом и его применение к задаче управления экономическим объектом // Вестник Томского государственного университета. Управление, вычислительная техника и информатика. 2012. № 1(18). С. 5-15.
9. Киселева М.Ю., Смагин В.И. Управление с прогнозирующей моделью с учетом запаздывания по управлению // Вестник Томского государственного университета. Управление, вычислительная техника и информатика. 2010. № 2(11). С. 5-12.
10. KiselevaM.Y., Smagin V.I. Model predictive control of discrete systems with state and input delays // Вестник Томского государственного университета. Управление, вычислительная техника и информатика. 2011. № 1(14). С. 5-12.
11. Kalman R.E., Busy R. A new results in linear filtering and prediction theory // Trans. ASME J. Basic Engr. 1961. V. 83. P. 95-108.
12. Амосов А.А., Колпаков В.В. Скалярно-матричное дифференцирование и его применение к конструктивным задачам теории связи // Проблемы передачи информации. 1972. № 1. С. 3-15.
Приступа Марина Юрьевна ООО «Битворкс» (г. Томск)
Смагин Валерий Иванович
Томский государственный университет
E-mail: kiselevamy@gmail.com; vsm@mail.tsu.ru Поступила в редакцию 25 февраля 2013 г.
Pristupa Marina Yu., Smagin Valery I. (Bitworks software company, Tomsk State University). Model predictive control of the output time-varying discrete systems with constraints on the control.
Keywords: discrete time-dependent systems, model predictive control, output control.
The problem of synthesis of predictive control, built on the basis of monitoring the output of the system for non-stationary linear object within the constraints on the control is considered. Forecasting is based on the calculation of time-dependent state estimates of the object using Kalman extrapolation and estimates of an unknown input.