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

Г.Н. Решетникова

СЛЕДЯЩАЯ СИСТЕМА АДАПТИВНОГО УПРАВЛЕНИЯ С ПРОГНОЗИРУЮЩЕЙ МОДЕЛЬЮ ПОНИЖЕННОГО ПОРЯДКА

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

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

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

1. ПОСТАНОВКА ЗАДАЧИ

Пусть математическая модель объекта управления задана в виде

X(t ) = A(t) x(t) + B (t)u(t) + F (t )q(t ),

x(to) = Xo,

U(t) = u(t), u(t0) = u0 , (1)

где x(t) є Rn - вектор состояния; u(t), u(t) є Rm -вектор положения органов управления и вектор управления, характеризующий скорость изменения u(t); A(t ), B (t), F (t) - матрицы модели объекта; q(t) -вектор внешних возмущений, которые описываются белым гауссовским шумом с характеристиками M {q(t )} = q (t ),

M {(q(t ) - q (t ))(q(T) - q (t))t } = Q(t )S(t - t) .

Пусть информация о состоянии объекта поступает в дискретные моменты времени к, соответствующие tk = t0 + кAt, к = 0,1,..., а управления являются кусочно-постоянными функциями на каждом интервале выдачи управляющих воздействий:

u(t) = и(кX tk <t < tk+1, tk+1 - tk =At.

Будем предполагать, что математическая модель информационной системы имеет вид

y (k) = Hx(k) + r (к), (2)

где y(k) є R1 - вектор, содержащий информацию о состоянии объекта; H - матрица, характеризующая

полноту информации об объекте; г (к) - погрешности информационной системы, которые будем считать гауссовскими последовательностями шумов с характеристиками:

М{г(к)} = 0, М{г(к)гТ О)} = Я8к,] , х(к) - вектор состояния дискретной модели, эквивалентной (1).

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

х(к +1) = А (к, 9(к)) х(к) +

+В (к )и(к) + Г (к )q(k), х(0) = х0,

и (к +1) = и(к) + Д/и(к), и (0) = и0, (3)

где 9(к) є Ям - вектор неизвестных параметров модели объекта, при этом предполагается, что элементы матриц А(-),В(-) - линейно зависят от компонент вектора 9(), т.е.

А (к) = 1п +д£А((к),

В(к) = Д/В (£к),

Г (к) = ТДг (/к).

Здесь 1п - единичная матрица п-го порядка; q(k) - гауссовская последовательность шумов с характеристиками

М {ц(к)} = Ц (к),

М{(ц(к) -Ц(к))(ц(]) -Ц(]))Т } = Q(k)Ь] .

Кроме того, будем предполагать, что априорные распределения векторов х0 и 90 являются гауссовскими:

М {х0} = х0,

М{(х0 - х0 )(х0 - х0)Т } = Рх0 ,

М {90} = 00,

М{(00 -00)(00 -00)Т} = Р90 и М{г(к)цт (])} = 0.

Формирование следящей системы будем осуществлять на основе минимизации критерия обобщенной работы [1]:

Г/к + ^Д/

3 = 0,5М | | [(х(/) - хг (/к ))Т С(х(/) - хг (/к)) +

I ‘к

+иТ (/)Ди(/) + и^р(/) + иТ (/)^(ОМ, (4)

где С, Б2 - неотрицательно определенные, а Д - положительно определенная весовые матрицы; 1р Д/ -

длина скользящего интервала оптимизации; х (•) -вектор заданных состояний, за которым осуществляется слежение; иор (•) - оптимальное управление, доставляющее минимум функционалу (4). При этом поведение объекта на интервале оптимизации [(к, (к + /рД(] будем описывать прогнозирующей моделью:

*м (У +1) = А(к, 0(к)) *М (У) +

+ В(к, 0(к ))и (к) + Е (к )д (к), хм (У = к) = х(к), у = к, к +1,..., 1р -1, (5)

где м - указывает на принадлежность прогнозирующей модели; х()0() - оценки векторов состояния и параметров, определяемые с помощью фильтров Кал-мана [2].

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

и(к) = и0р (к) = Д-'И^к), (6)

где И2(к) - да-мерный вектор, который является решением в обратном времени системы обыкновенных дифференциальных уравнений:

И (к + 1р - (у +1)) = Ат (к, 0(к))И1 (к + - У) +

+Д((ом (к + 1р - у) - хг ((к)), И1(к + 1р) = 0,

^(к + 1р - (у +1) = ^(к + /р - у) +

+Вт (к, 0(к))И1 (к + /р - у) + Д(Д2и(к), И2(к + /р) = 0, ом (к + /р - (У +1)) = 2ом (к + /р - у) -

- А (к, 0(к ))ом (к + /р - у) - В(к, 0(к ))и(к) - Е (к )д (к),

ом(к + /р) = хм(к + /р),у = к,к +1,...,к + /р -1. (7)

Последнее уравнение в (7) описывает прогнозирующую модель на интервале оптимизации в обратном времени.

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

х(() = Ах(() + Ви ((), х((0) = х0 требуется найти параметры агрегированной системы хр (() = Архр (() + Ври((), хр ((0) = х0р порядка р < п , для которой выполняется соответствие:

хр (() = НрХ((),

где Нр - матрица агрегирования размерности р х п ранга р. При этом последнее соотношение выполнено, если

Ар (() = НрА((), Вр = НрВ, Хор = НрХо .

2. АГРЕГИРОВАНИЕ ЛИНЕЙНОЙ ДИСКРЕТНОЙ СТОХАСТИЧЕСКОЙ СИСТЕМЫ

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

Для агрегирования воспользуемся дискретной системой, определенной на интервале [(к, (к + /р М]:

х( у +1) = А (к) х( у ) + В (к )и (к) + Е (к )д( у), х(] = к) = х(к),у = к,к +1,...,к + /р -1, (8)

где д( у) - гауссовская последовательность шумов с характеристиками:

М {д( У)} = д (к),

М{(д(У) - д(У))(д(/) - д(1))т } = 0(к)8у,.

Будем предполагать, что

М{х(у)} = х(У),У = к,к +1,...,к + /р -1,

х (У +1) = А(к) х (у) + В(к )и (к) + Е (к )д (к), х(у = к) = х(к),у = к,к +1,...,к + /р -1. (9)

Систему пониженного порядка р < п с вектором состояния хр (•) для (8) будем строить в виде

хр (У +1) = Ар (У) хр (У) + НрВ(к )и (к) + НрЕ (к )д( у ),

хр(у = к) = Нрх(к),у = к,к +1,...,к + /р -1, (10)

где хр (•) е ,Кр - вектор состояния системы пониженного порядка; Н р - матрица агрегирования размерности р х п ранга р;

Ар (у) = НрА(к) Д( у).

Здесь Д( у) - матрица оператора обратного преобразования Д : Яр ^ Яп размерности п х р , которую будем строить таким образом, чтобы восстановление состояния исходной системы по состоянию агрегированной

х() = Д{хр (•)}

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

Д( У) = Н++ Нр^ (У), (11)

где Н+ - матрица, псевдообратная к Н ;

Н = 1п - Н++ Н - проекционная матрица; 5^(у) - произвольная матрица размерности п х р . Так как матрица Д( у) вида (11) является решением уравнения НрД(У) = 1р , что проверяется непосредственно подстановкой, то НрД(У)хр(У) = хр(У), У = к,к +1,...,к + /р.

Матрицу £(у) будем определять из условия минимума функционала:

./(к) = М{(х(у) - х(у ))т (х(у) - х(у))} =

= 1гМ{((х(у) - Д( у) хр (у))(х(у) - Д( у) хр (у))т }. (12)

Ошибка аппроксимации системы (8) системой (10) определяется соотношениями

е( У +1) = х( У +1) - х( У +1) =

= Д(у + 1)НрА(к)е(у) + [ 1п - Д(У +1)Нр]х(у +1), (13) е(у = к) = [ 1п - Д(У = к)Нр ]х(к), (14)

а уравнение для ковариационной матрицы ошибки имеет вид

Ре (У +1) =

= Д(У +1)Нр!(У +1)Нтр Дт (У + 1) - Д(У +1)Нр¥т (У + 1) -

-V(У +1)Нтр Дт (У +1) + Ох (У +1), (15)

где Ц у +1) = Ох (у +1) + А(к) Ре (У) Ат (к) -

- А(к )Оех (У) Ат (к) - А (к О (у )Ат (к); (16)

V (у +1) = Ох (у +1) - А (к О (у )Ат (к), (17)

( (•)= М {х()хт (•)}).

В предположении, что векторы х(у) и д(у) не коррелированные, получим

Ох (у +1) = А (к)Ох (у) Ат (к) +

+А (к) х( у)[ В (к )и( к) + Е (к )д(к )]т +

+[ В(к )и( к) + Е (к )д(к )]хт (у) Ат (к) + С (к); (18)

С (к) = Е (к )0(к) Ет (к) +

+[В(к)и(к) + Е(к)д(к)][В(к)и(к) + Е(к)д(к)]т . (19) Величина Оех (•) при некоррелированных е(у) и д(у) определяется следующим образом:

Оех (У +1) = М {е( У +1) х( (У +1)} =

= Ох(у +1) - Д(у +1)НрУт (у +1). (20)

Начальные условия для (15) - (20) имеют вид

Ох (У = к) = М{х(к)х (к)} = Охк ;

(21)

Оех (У = к) = [1п - Д(У = к)Нр ]Ох (У = к); (22)

Ре (У = к) = [ 1п - Д( У = к) Нр ] х хОх(у = к)[/п - Д(у = к)Нр]т . (23)

Подставляя в (15) выражение (11) для Д(у +1) и приравнивая градиент 1г{Ре (у +1)} по £(у +1) к нулю, получим уравнение для определения £(у +1):

НрБ(У +1)НрЬ(У +1)Нр = НрУ(У +1)Нтр . (24)

Если гапк{Цу)}> р для у = к,к +1,...,к + /р, то

получаем уравнения, единственным образом определяющие матрицу Д( у +1) в виде

Д(У +1) = Н+р + HpV(У +1)Нтр [Нр!(У +1)Нтр ]-1; (25) Д( У = к) = Ох (У = к )Нр [НрОх (У = к )Нр ]-1. (26)

При этом, так как для оптимальной матрицы Д( у), определяемой согласно (25) - (26) выполняется соотношение Ре(у) = Оех(у), у = к,к +1,...,к + /р, что достаточно легко доказывается по индукции, то

Д(У +1) = Цу +1)Нтр [Нр!(у +1)Нтр ]-1 ; (27)

Ц у +1) = Ох (у +1) - А (к) Ре (у) Ат (к); (28)

Ре (у +1) = Ох (у +1) - Д(у + 1)НрЦу +1). (29)

Таким образом, построение системы пониженного порядка на интервале [(к, (к + /р Д( ] сводится к вычислению установившегося значения матрицы Ар (•) системы (10) по рекуррентным формулам:

Ох (У +1) = А (к )Ох (У) Ат (к) +

+А( к) х( у)[ В(к )и (к) + Е (к )д(к )]т +

+[ В(к )и (к) + Е (к)д (к )]хт (у )Ат (к) + С (к),

Ц У +1) = Ох (У +1) - А (к) Ре (У) Ат (к),

Д(У +1) = Ду +1)Нтр [Нр!(у +1)Нтр ]-1,

Ар (у +1) = НрА(к) Д( у +1),

Ре (у +1) = Ох (у +1) - Д( у +1) Нр!( у +1), х (у +1) = А (к) х (у) + В(к )и (к) + Е (к )д(к), (30)

(у = к,к +1,...,к + /р -1), при следующих начальных и постоянных значениях:

Ох (У = к) = Охк,

Д(у = к) = Ох (у = к)Нтр [НрОх (У = к)Нтр ]-1,

Ре (У = к)=[ 1п - Д (У = к) Нр ]Ох (у = к)[ 1п - Д (у = к) Нр ]т, х (у = к) = х(к),

С (к) = Е (к )0( к) Ет (к) +

+[В(к)и(к) + Е(к)д(к)][В(к)и(к) + Е(к)д(к)]т . (31)

Установившееся значение Ар (•), обозначаемое Ар (к), определяется при выполнении условия

Ар (к) =

Ар (У), если

(Г {Ре (У)}- (Г {Ре (У - 1)}

]</р

(32)

(г {Ре (У -1)}

Акр (/р), в противном случае,

и система пониженного порядка в соответствии с (10) будет иметь вид

хр (У +1) = Ар (к)хр (у) + Вр (к)и(к) + Ер (к)д(к),

хр (у = к) = Нрх(к),у = к,к +1,...,к + /р -1, (33)

где е - заданная точность агрегирования,

Вр (к) = НрВ(к), Ер (к) = НрЕ(к).

3. СИНТЕЗ АДАПТИВНОГО УПРАВЛЕНИЯ С ПРОГНОЗИРУЮЩЕЙ МОДЕЛЬЮ ПОНИЖЕННОГО ПОРЯДКА

Будем формировать управляющие воздействия для (1) на основе минимизации функционала обобщенной работы (4) при описании поведения объекта на интервале оптимизации [(к, (к + /рД(] прогнозирующей моделью пониженного порядка [5]:

хМ (У+1)=Ар (к,0(к))хМ (у )+Вр (к,0(к))и(к)+Ер (к)д(к), хр(У = к) = Нрх(к), у = к,к +1,...,к + /р -1, (34)

где матрица динамики Ар (к, 0(к)) определяется согласно (30) - (33) при замене А(к) на А(к,0(к)), В(к) на В(к, 0(к)), х(к) на х(к). Тогда, если в (4) вместо х() воспользоваться восстановленным состоянием х(•), а затем обозначить

хр (к) = Нрхг (к), Ср = Д(к)СДт (к), где Д(к) - матрица оператора обратного преобразования, определяемая одновременно с Ар (к, 0(к)), то управление будет иметь вид

и(к) = - Д-1И2р (к). (35)

В (35) да-мерный вектор Ир (к) является решением в обратном времени системы обыкновенных дифференциальных уравнений:

Ир (к + /р - (у +1)) = А'т (к, 0(к))И1р (к + /р - у) +

+МК (к + /р - у) - хгр ((к)), И1р (к + /р) = 0,

Ир (к + /р - (У +1) = Шр (к + /р - У) +

+Вт (к, 0(к))И1р (к + /р - у) + Д(Д2и(к), И2р (к + /р) = 0, ом (к + /р - (у +1)) = 2ом (к + /р -1) -

- Ар (к ,0(к ))ом (к+/р - у) - Вр (к ,0(к ))и (к) - Ер (к )д (к),

°м (к + /р ) = хмм (к + /р ), У = к, к + 1,..., к + /р -1 . (36)

4. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ

Численное моделирование осуществлялось для экономической модели [6], описывающей процесс производства, хранения и сбыта продукции, в виде

системы обыкновенных дифференциальных уравнений (1). При этом компоненты вектора состояния х(() = (х1 ((), х2 ((), х3 (())т имеют следующий смысл: х1 ((), х2 (() - объемы продукции на рынке и у потребителя; х3(() - прибыль от реализации; и(() - темп производства; и(() - управление, которое задает изменение темпа производства; д(() - вектор, описывающий действия случайных факторов, который будем считать нормальным гауссовским шумом. Матрицы А(() и В (() определяются следующим образом:

( -пс (ф7 (() -^ (()) -псфх (() 0^

пс (ф7 (() -Фv (()) -псФх (() - к 0

V спс (ф7 (() -Фv (()) - к2 спс фх (() 0,

1

A(t) =

5 (t) = 11 0 —

где к1, к2, с, пс - коэффициенты, характеризующие темп потребления, плату за хранение, торговую накрутку, скорость продажи; ф7 ((), фх ((), фV (() - функции, описывающие математические модели потенциального спроса, ситуацию на рынке и поведение покупателя.

Моделирование осуществлялось по данным балансовых счетов предприятия, производящего продукты питания при наличии неизвестных параметров, образующих вектор 0() = (к1, фг (•), фх (•), фV (•))т . Синтез управлений осуществлялся при слежении за состоянием, которое соответствует максимально возможной прибыли для рассматриваемого предприятия при существующей ситуации на рынке. Значение максимально возможной прибыли определялось с помощью элементов факторного анализа и аппарата производственных функций типа Кобба - Дугласа [7].

Управление формировалось по информации, характеризующей только динамику изменения прибыли (Н = (0 0 1)), а агрегирование прогнозирующей модели осуществлялось с помощью матрицы Нр = (0 0 1), что соответствует слежению за величиной хр (•), которая в данном случае совпадает со значением максимально возможной прибыли. В результате моделирования получилось, что наращенная прибыль практически совпадает с отслеживаемой.

ЛИТЕРАТУРА

1. Решетникова Г.Н. Адаптивное управление линейными дискретными стохастическими системами по критерию обобщенной работы // Изв. АН СССР, Технич. кибернет. 1989. № 3. С. 196.

2. Брамер К., ЗиффлингГ. Фильтр Калмана - Бьюси. М.: Наука, 1982.

3. AokiM. Control of large-scale dinamie sistems by aggregation // IEEE Trans. Automat. Control. 1968. V. AC-13. No. 3. P. 246 - 252.

4. Домбровский В.В. О приближенном агрегировании линейных стохастических систем // Автоматика и телемеханика. 1989. № 7. С. 102 -109.

5. Решетникова Г.Н. Построение и использование прогнозирующих моделей пониженного порядка для синтеза цифрового адаптивного управления. Деп. В ВИНИТИ № 1556 - В96 от 15.05.96.

6. Решетникова Г.Н. Синтез следящей системы адаптивного управления темпом производства по локальному критерию // Вестник ТГУ. 2004. № 284. С. 166 - 168.

7. Терехов Л.Л. Производственные функции. М.: Статистика, 1974.

Статья представлена кафедрой прикладной математики факультета прикладной математики и кибернетики Томского государственного университета, поступила в научную редакцию «Кибернетика» 1 июня 2005 г.