Математические модели механических систем. Построение и моделирование линейных математических моделей систем управления

1 Курс «Математическое моделирование технических объектов и систем управления»-2016. Лекция 2. Математические модели механических систем. Построение и моделирование линейных математических моделей систем управления Математические модели механических систем. Линеаризованные математические модели первого приближения. Математический аппарат анализа линейных моделей. Методы синтеза линейных моделей. Краткое описание функций пакета Control System Toolbox. Наблюдатели полного и сокращенного порядков. Математические модели механических систем. Дифференциальные уравнения движения систем материальных точек. Рассмотрим движение системы материальных точек Pk , k  1, 2,..., N относительно некоторой инерциальной системы координат. Будем считать, что на положения и скорости точек системы наложены ограничения геометрического или кинематического характера, называемые связями. Системы, с такого рода связями, называются несвободными. Аналитически связь для каждой точки системы обычно выражается уравнением вида: f (t , rk , rk )  0, k  1, 2,..., N , где rk - радиус вектор, характеризующий положение точки Pk , а rk  vk - скорость данной точки. Если скорость rk не входит в уравнение, то есть связь имеет вид: f (t , rk )  0 , то такую связь называют геометрической. В общем случае, связь называется дифференциальной или кинематической. Примем, что каждая кинематическая связь может описываться линейными уравнениями вида: N l k 1 k  rk  D  0 , где векторы lk и скаляр D представляют собой заданные функции от параметров t , rk . Будем предполагать, что каждая геометрическая связь f (t , rk )  0 влечет за собой, как следствие, дифференциальную связь вида: N f  r k 1 k  rk  f  0, t которая будет эквивалентна конечной связи f (t , rk )  ck . Такие связи будем называть интегрируемыми. В прямоугольных декартовых координатах уравнение связи для каждой точки можно записать в форме: f (t , xk , yk , zk , xk , yk , zk )  0 f (t , xk , yk , zk )  0 Или, для линейных связей, в форме: N (A k 1 N  xk  Bk  yk  Ck  zk )  D  0 k f  ( x k 1 k  xk  f f f  yk   zk )  0 yk zk t 2 Такие связи называются стационарными, если выполняются условия D  0 , f 0 и t функции Ak , Bk , Ck не зависят явно от параметра t . Система материальных точек называется голономной, если она является свободной (то есть не имеющей связей) или имеет геометрические и/или дифференциальные интегрируемые связи. При наличии дифференциальных неинтегрируемых связей система называется неголономной. Система называется склерономной, если на нее наложены только стационарные связи. В противном случае она называется реономной. Пример. Пусть две материальные точки постоянной длины P1  ( x1 , y1 , z1 )T , P2  ( x2 , y2 , z2 )T соединены стержнем l . Тогда уравнения связи имеет вид: ( x1  x2 ) 2  ( y1  y2 ) 2  ( z1  z2 ) 2  l 2  0 . Таким образом, такая система является голономной склерономной системой. Если точки соединены стержнем переменной длины l (t ) , то уравнение связи примет форму: ( x1  x2 ) 2  ( y1  y2 ) 2  ( z1  z2 ) 2  l (t ) 2  0 . Это будет голономная реономная система. Пусть теперь две материальные точки могут двигаться только по плоскости, и соединены стержнем постоянной длины l . При этом точки могут двигаться только так, чтобы скорость середины стержня была направлена вдоль стержня. Тогда уравнения связей имеют форму: z1  z2  0 , ( x1  x2 ) 2  ( y1  y2 ) 2  l 2  0 , x1  x2 y1  y2  x1  x2 y1  y2 . Тогда такая система будет неголономной, так как последнее уравнение определяет дифференциальную неинтегрируемую связь. Возможные и виртуальные перемещения. Пусть на материальную систему наложены d геометрических связей fi (t , rk )  0 , i  1, 2,..., d ; k  1, 2,..., N и s дифференциальных связей N l j 1 r, j  rj  Dr  0 , r  1, 2,..., s . Заменим геометрические связи вытекающими дифференциальными связями получим систему: N l j 1 N r, j fi  r k 1 k  v j  Dr  0 , r  1, 2,..., s  vk  fi  0 , i  1, 2,..., d , t где vk  rk . Тогда систему векторов vk будем называть возможными скоростями для некоторого момента времени t , если она удовлетворяет d  s дифференциальным уравнениям. То есть возможные скорости – это скорости, допускаемые связями. Очевидно, 3 что в общем случае, в момент времени t существует бесчисленное множество систем возможных скоростей, хотя при действительном движении системы реализуется только одна из этих систем скоростей. Определим систему возможных бесконечно малых перемещений (или просто систему возможных перемещений) как: drk  vk  dt , k  1, 2,..., N . Такие перемещения будут удовлетворять уравнениям вида: N l j 1 N r, j  drj  Dr  dt  0 , r  1, 2,..., s fi  r k 1  drk  k fi  dt  0 , i  1, 2,..., d t Рассмотрим две системы возможных перемещений drk  vk  dt , drk  vk  dt , k  1, 2,..., N . Назовем разность rk  rk  rk виртуальным перемещением. Очевидно, что система виртуальных перемещений будет удовлетворять однородным уравнениям: N l j 1 N r, j fi  r k 1  rj  0 , r  1, 2,..., s  rk  0 , i  1, 2,..., d k Таким образом, виртуальные перемещения совпадают с возможными перемещениями при «замороженных» связях. То есть, когда время t фиксируется и связи «застывают» в той конфигурации, которую они имели к данному моменту времени ( dt  0 ). Можно сказать, что виртуальные перемещения представляют собой перемещение возможного положения системы точек в момент t в другое бесконечно близкое возможное положение. Очевидно, что при стационарных связях виртуальные перемещения совпадают с возможными перемещениями. В декартовых координатах уравнения, определяющие виртуальные перемещения могут быть записаны в форме: N (A j 1 N r, j fi  ( x k 1 k  x j  Br , j  y j  Cr , j  z j )  0, r  1, 2,..., s  xk  fi f  yk  i  zk )  0, i  1, 2,..., d yk zk Если эти d  s уравнений независимы, то среди виртуальных приращений будет n  3N  d  s независимых величин. Число n называется числом степеней свободы системы материальных точек. Предположим, что в точках Pk , k  1, 2,..., N системы, имеющих массу mk , k  1, 2,..., N , приложены активные силы Fk , k  1, 2,..., N . Если бы связи отсутствовали, то имело бы место соотношение: mk  wk  Fk , k  1, 2,..., N , 4 где wk - соответствующее ускорение точки Pk . Однако, при наличии связей, величины wk могут быть несовместимы с соответствующими дифференциальными уравнениями связей. Действительно, продифференцировав уравнения по времени получим: N N dlr , j j 1 j 1 dt  lr , j  w j    vj  dDr  0 , r  1, 2,..., s dt N fi d fi d fi  wk   ( )  vk   0 , i  1, 2,..., d ,  dt t k 1 rk k 1 dt rk N В общем случае, система ускорений wk , k  1, 2,..., N может не удовлетворять этим уравнениям. Таким образом, связи будут создавать некоторые дополнительные силы Rk , k  1, 2,..., N , которые будут действовать на точки системы. Эти силы называются реакциями связей. При этом реальные ускорения точек системы будут удовлетворять уравнениям: mk  wk  Fk  Rk , k  1, 2,..., N . Активные силы Fk обычно задаются как некоторые известные функции вида: Fk  Fk (t , r , v), r  (r1 , r2 ,..., rN )T , v  (v1 , v2 ,..., vN )T , k  1, 2,..., N Тогда основная задача исследования динамики несвободной системы имеет следующую форму – по заданным активным функциям Fk  Fk (t , r , v ), k  1, 2,..., N и совместимым со связями начальному положению r и начальным вектором скорости v определить движение системы и реакции связей Rk , k  1, 2,..., N . В координатной форме задача сводится к нахождению 6N величин ( xk , yk , zk , Rx , k , Ry ,k , Rz ,k ) . При этом число уравнений будет равно 3N  d  s  6N , то есть система имеет 6N  (3N  d  s)  3N  d  s  n независимых соотношений между величинами или степеней свободы. Чтобы определить эти независимые величины обычно ограничиваются классом идеальных связей, который задается условием: N R k 1 k  rk  0 , или в развернутом виде, как: N  (R k 1 x,k  xk  Ry , k  yk  Rz , k  zk )  0 . То есть принимается, что сумма работ реакций этих связей на любых виртуальных перемещениях всегда равна нулю. С помощью данного соотношения можно выразить 3N  n зависимых приращений (xk , yk , zk ) через n независимых приращений и приравнять нулю коэффициенты при этих независимых приращениях. Таким образом, основная задача динамики становится определенной. 5 Пример. Пусть материальная точка вынуждена двигаться по подвижной гладкой поверхности, как это показано на рисунке 2.1. Рис.2.1 В этом случае возможная скорость материальной точки, а значит и возможное малое перемещение dr  v  dt будет лежать не в касательной плоскости. Однако, виртуальное перемещение r , которое представляет собой бесконечно малое возможное перемещение для «замороженной» поверхности будет лежать в касательной плоскости. В этом случае реакция «замороженной» подвижной гладкой поверхности будет направлена по нормали к поверхности, то есть будет выполняться соотношение как R  dr  0 ). То есть данная связь является идеальной. R  r  0 (в то время Общее уравнение динамики. Уравнения Лагранжа первого рода. Уже отмечалось, что материальные точки системы удовлетворяют уравнению mk  wk  Fk  Rk , k  1, 2,..., N . Если связи идеальные, то выполняется соотношение N R k 1 k  rk  0 . Отсюда получим следующее общее уравнение динамики несвободной системы под действием активных сил: N  (F k k 1  mk  wk )  rk  0 . Данное уравнение выражает необходимое и достаточное условие того, что движение, совместимое со связями, соответствует системе активных сил Fk  Fk (t , r , v), k  1, 2,..., N . Найдем выражение для реакций Rk , k  1, 2,..., N с помощью метода неопределенных множителей, предложенным Лагранжем. Для виртуальных перемещений точек системы будут справедливы соотношения: N l j 1 N r, j fi  r k 1  rj  0 , r  1, 2,..., s  rk  0 , i  1, 2,..., d k Умножая каждое равенство на произвольные скалярные множители  i , r и, складывая данные соотношения, получим: 6 N d k 1 i 1  ( Rk   i  s fi    r  lr , k )  rk  0 . rk r 1 Так как величина rk является произвольной, то получим: d Rk    i  i 1 s fi    r  lr , k . rk r 1 Подставляя полученное выражение в общее уравнение динамики несвободной системы, и, учитывая уравнения связей, найдем: d mk  wk  Fk    i  i 1 s fi    r  lr , k , k  1, 2,..., N rk r 1 fi (t , rk )  0 , N l j 1 r, j  v j  Dr  0 , где i  1, 2,..., d , r  1, 2,..., s . Данная система называется уравнениями Лагранжа первого рода. Таким образом, получим систему из 3N  d  s скалярных уравнений с 3N  d  s координатными неизвестными скалярными величинами ( xk , yk , zk , i , r ) , что позволяет найти величины реакций и уравнения движения. Однако, интегрировать такую систему уравнений крайне затруднительно из-за большого числе уравнений. Поэтому данные уравнения достаточно редко используются в приложениях. Принцип виртуальных перемещений Даламбера. Будем называть положением равновесия такое положение системы, в котором система будет находиться все время, если в начальный момент времени она находилась в этом положении и скорости всех ее точек равны нулю. Очевидно, для этого положения системы будет выполняться соотношение: N F k 1 k  rk  0 . Данное равенство представляет собой принцип виртуальных перемещений. То есть, для того, чтобы некоторое, совместимое со связями, положение системы было положением равновесия, необходимо и достаточно, чтобы в этом положении сумма работ активных сил на любых виртуальных перемещениях системы равнялась нулю. Чаще всего принцип виртуальных перемещений используют для систем со стационарными связями. Очевидно, что указанное уравнение представляет собой частный случай общего уравнения динамики, рассмотренное в окрестности положения равновесия, если к активным силам Fk  Fk (t , r , v), k  1, 2,..., N причислить фиктивные силы инерции  mk  wk , k  1, 2,..., N . Таким образом, можно сформулировать принцип Даламбера. 7 При движении системы любое ее положение можно рассматривать как положение равновесия, если к активным силам, действующим на систему в этом положении, прибавить фиктивные силы инерции. Данный принцип позволяет перенести приемы и методы решения статических задач на задачи динамики. Действительно, в положении равновесия реакции связей определяются соотношением: Fk  mk  wk   Rk , k  1, 2,..., N . То есть динамическая задача заменяется статической задачей, решения которой совпадают с исходной динамической задачей. Пример. Горизонтальный однородный вал равномерно вращается с угловой скоростью  . Перпендикулярно к оси вала на равных расстояниях от подшипников на вал эксцентрично насажен однородный диск. Найти давление на подшипники при вращении вала, имеющего массу M a . Схема такого устройства приведена на рисунке 2.2. Рис. 2.2 Равнодействующая сил инерции, действующих на отдельные элементы определяются соотношением: R  2   r  dm  M d  2  rc где dm , равные 2  r  dm , , Md - масса диска, rc - расстояние между геометрическим центром диска и осью вала. На основе принципа Даламбера определим статические давления на подшипники. На ось вала приложены следующие силы: сила веса вала инерционная сила M a  g , сила веса диска M d  g и R  M d  2  rc . Тогда давление на каждый подшипник определяется соотношением: N 1 1  (M a  M d )  g   M d  2  rc . 2 2 Сила N имеет максимальную величину в том положении, когда геометрический центр диска будет расположен под точкой О. Независимые координаты и обобщенные силы. Рассмотрим голономную систему из N материальных точек Pk , k  1, 2,..., N с геометрическими связями fi (t , r )  0 , i  1, 2,..., d . Из уравнений связи можно выразить d координат, как функции 3N  d остальных и времени t . Однако можно все 3N декартовых 8 координат выразить в виде функций от n  3N  d независимых параметров q1 , q2 ,..., qn и от параметра t : xk  k (t , q), yk   k (t , q), zk   k (t , q), q  (q1 , q2 ,..., qn )T . Данные равенства эквивалентны векторным равенствам rk  rk (t , q ), k  1, 2,..., N . Величины q1 , q2 ,..., qn называются независимыми обобщенными координатами. Каждой координате q j соответствует своя обобщенная сила Q j , j  1, 2,..., n . Для определения этих сил рассмотрим элементарную работу активных сил на виртуальных перемещениях: N A   Fk  rk , k 1 где rk  n rk  q j 1  q j , k  1, 2,..., N . Отсюда найдем, что элементарная работа может быть j записана с помощью соотношения: n N n rk rk A   Fk   q j   ( Fk  )  q j   Q j  q j , q j k 1 j 1 q j j 1 k 1 j 1 N n где обобщенные силы определяются равенствами Q j  N F k 1 k  rk , j  1, 2,..., n . Но q j приращение q j независимых координат q j могут быть произвольными. Тогда можно записать систему равенств: Q j  0, j  1, 2,..., n . Отсюда получим, что положение голономной системы является положением равновесия только в том случае, когда в этом положении все обобщенные силы равны нулю. Пример. Рассмотрим твердое тело, которое может только вращаться вокруг некоторой неподвижной оси. A  M   , где M - суммарный момент всех активных сил относительно оси вращения. Отсюда найдем, что Q  M . Соответствующий угол  может быть взят в качестве независимой координаты. Тогда Уравнения Лагранжа второго рода. Рассмотрим общее уравнение динамики свободной голономной системы: N  (F k 1 k  mk  wk )  rk  0 . Элементарная работа активных сил в такой системе в независимых координатах имеет вид: N n k 1 j 1 A   Fk  rk   Q j  q j , 9 где Q j  N F k 1 k  rk , j  1, 2,..., n . q j Элементарная работа сил инерции  mk  wk , k  1, 2,..., N может быть записана в форме: N n k 1 j 1 Ai   mk  wk  rk   Z j  q j , где Z j  N  mk  wk  k 1 N rk dr r   mk  k  k . q j k 1 dt q j После ряда преобразований обобщенные силы инерции можно преобразовать к следующему виду: N rk r d N d T T Z j   mk  rk    mk  rk  k    , j  1, 2,..., n dt k 1 q j k 1 q j dt q j q j 1 N 2 где T   mk  rk - кинетическая энергия системы. 2 k 1 Из общего уравнения динамики получим, что A  Ai  0 или N  (Q j 1 j  Z j )  q j  0 . Отсюда получим, что уравнения динамики системы могут быть записаны в виде: d T T    Q j , j  1, 2,..., n . dt q j q j Данная форма уравнений называется уравнениями Лагранжа в независимых координатах или уравнениями Лагранжа второго рода. Величины q j называются обобщенными скоростями. В случае несвободной системы необходимо также определить реакции связей Rk , k  1, 2,..., N , которые не входят в уравнения Лагранжа. После того, как уравнения Лагранжа будут проинтегрированы и найдены функции q j (t ), j  1, 2,..., n можно найти функции rk (t ), vk  rk (t ), wk  rk (t ), Fk (t , r , r ), r  (r1 , r2 ,..., rN ) T . Неизвестные реакции находятся из соотношения: Rk  mk  wk  Fk , k  1, 2,..., N . Пример. Рассмотрим твердое тело, которое может только вращаться вокруг некоторой неподвижной оси. Соответствующий угол QM , где M  может быть взят в качестве независимой координаты. Отсюда найдем, что - вращающий момент. С другой стороны T 1  I  2 , где I 2 - момент инерции тела 10 относительно оси вращения. Уравнения Лагранжа после подстановки принимает вид T T  I ,  0, Q  M   I   M . Канонические уравнения Гамильтона. Пусть обобщенные силы Q j  0, j  1, 2,..., n являются потенциальными, то есть существует потенциал сил (потенциальная энергия)   (t , q) такой, что Q j   Тогда уравнения Лагранжа  . q j d T T     , j  1, 2,..., n можно записать в форме: dt q j q j q j d L L    0, j  1, 2,..., n , dt q j q j где L  T   . Функция L называется функцией Лагранжа или кинетическим потенциалом. Кинетический потенциал, также же как и кинетическая энергия T , представляет собой функцию второй степени относительно обобщенных скоростей: L  L2  L1  L0 , n 1 n где L2   ci , j  qi  q j , L1   ci  qi , L0  c0 . Коэффициенты ci , j , ci , c0 являются 2 i , j 1 i 1 функциями от координат q1 , q2 ,..., qn и времени t . Для систем общего вида, движение которых описывается уравнения Лагранжа, кинетический потенциал может описываться произвольной функцией 2 L L  L(t , q, q ), q, q   , для которой предполагается выполнение условия det( )  0. qi q j n Переменные t , q, q , через которые выражается функция L , часто называют переменными Лагранжа. Гамильтон предложил в качестве основных переменных, характеризующих состояния системы, выбрать величины t , q, p , где p  ( p1 , p2 ,..., pn ) . Здесь pi - обобщенные T импульсы, определяемые равенством: pi  L , i  1, 2,..., n . qi Переменные t , q, p называются переменными Гамильтона. Так как гессиан функции L отличен от нуля, то можно выразить обобщенные скорости через обобщенные координаты и обобщенные импульсы: qi   i (t , q, p ), i  1, 2,..., n . 11 Таким образом, любая функция от переменных Лагранжа может быть преобразована в функцию от переменных Гамильтона: G (t , q, q)  G (t , q, p) . Введем в рассмотрение функцию: n H (t , q, p)   pi  qi  L , i 1 где функции qi , L, i  1, 2,..., n зависят от переменных t , q, p . Гамильтон показал, что уравнения движения системы могут быть записаны в виде системы 2n дифференциальных уравнения первого порядка: dqi H dpi H  ,  , i  1, 2,..., n dt pi dt qi Данные уравнения называются каноническими уравнениями в форме Гамильтона. Если функция Гамильтона не зависит явно от параметра t , то есть: H  H (q, p) , то такая система называется консервативной. При этом H (q, p)  h  const , где величина h называется обобщенным интегралом энергии. Если система является склерономной, то можно показать, что функция Гамильтона равна H  T   , то есть представляет собой полную энергию, выраженную через переменные Гамильтона t , q, p . Пример. Горизонтальная рейка, вращается около вертикальной оси, а вдоль рейки движется груз массы Сила, действующая на груз, имеет потенциал  (r ) , где r m. - расстояние от центра масс груза до оси вращения. Необходимо найти уравнения движения и интеграл энергии.  I  m  d 2 момент инерции рейки относительно оси вращения. Введем следующие независимые координаты системы: , r . Тогда кинетическая энергия системы Обозначим через угол поворота рейки, а через может быть записана в виде: T 1 1 1  I  2   m  ( r 2  r 2   2 )   m  ( r 2  ( r 2  d 2 )   2 ) . 2 2 2 Отсюда найдем: pr  T  mr r , p  T  m  (r 2  d 2 )   .  Очевидно, что система является консервативной, так как полная энергия системы равна: p2 1 2 H T    ( pr  2 )   (r ) . 2m r  d2 Канонические уравнения Гамильтона имеют вид: 12 dr pr  , dt m p d  dt m  (r 2  d 2 ) r  p2 dpr    (r ), dt m  (r 2  d 2 )2 dp 0 dt Интеграл энергии равен: p2 1 2  ( pr  2 )  (r )  h 2m r  d2 Линеаризованные математические модели первого приближения. Линейные математические модели систем применяются в некоторой малой окрестности положения равновесия, когда управляющие сигналы ограничены по величине. В этом случае, исходная нелинейная модель может быть достаточно хорошо описана с помощью правильно построенной линейной модели. Если линейная модель будет построена плохо, то и результаты построения системы будут плохими, и использование такой системы управления с нелинейным объектом приведет к плохому качеству управления. Таким образом, процесс линеаризации нелинейных уравнений является очень важным при проектировании систем управления. Процедура линеаризации базируется на представлении нелинейных функций в виде ряда Тейлора, относительно некоторой точки равновесия. В дальнейшем, принимая малыми возможные отклонения переменных, мы ограничиваемся только линейными членами в ряде Тейлора, что позволяет построить линейную модель объекта. Линейная аппроксимация нелинейных математических моделей. Рассмотрим нелинейное уравнение вида: x(t )  F ( x(t ), u(t )) , где функция F определяет отображение F :     . Точка x  называется изолированной или гиперболической точкой равновесия, если существует такой постоянный n m n n   m , что F ( x , u )  0 , и в ближайшей окрестности точки ( x , u ) нет других точек F x  x равновесия, а det( |u u )  0 . Если, начальное состояние системы равно x(t0 )  x при входном x воздействии u , то, очевидно, что x(t )  x для t  t0 . входной сигнал u Определим вариации переменных относительно точки равновесия: x(t )  x(t )  x , u(t )  u(t )  u , x (t ) решение уравнения x(t )  F ( x(t ), u(t )) . Тогда можно записать следующее уравнение: x(t )  F ( x  x(t ), u  u(t )) . Применим разложение в ряд Тейлора для правой части уравнения где и, оставим только линейные его члены. Тогда можно записать приближенное равенство: x(t )  F ( x ), u )  Так как F ( x , u )  0 , то получим: F x  x F x  x |u u x(t )  |u u u x u 13 F x  x F x  x |u u x(t )  |u u u . x u F x  x B |u u nm являются постоянными, вещественными u x(t )  Матрицы A  F x  x |u u nn x и матрицами, которые также называются матрицами Якоби. Линейная стационарная система x(t )  A  x(t )  B  u называется линейной системой первого приближения. Используя такую x(t ) ее переменных и применяя малые управления u (t ) относительно точки равновесия ( x , u ) . модель можно спроектировать систему управления в области вариации Пример.  0 ): x  x(6  2 x  y) y  y(4  x  y) Здесь имеется две изолированные точки равновесия x1  y1  0 и x2  y2  2 . В окрестности первой точки x1  y1  0 матрица Якоби равна:  x  x 0  6 0   6  4x  y  | y  0    . J ( x1 , y1 )   4  x  2 y   y  0 4 Очевидно, что Якобиан: det( J ( x1 , y1 ))  0 . Тогда можно записать следующую систему первого  x   6 0  x    . приближения:      y   0 4  y  В окрестности второй точки равновесия x2  y2  2 матрица Якоби равна: 6  4x  y  x  x2   4  2 J ( x2 , y2 )   |  ; det( J ( x2 , y2 ))  0 .  y 4  x  2 y  y  2   2  2    x    4  2  x    . Тогда система уравнений первого приближения имеет вид:      y    2  2  y  Пусть нелинейная система имеет уравнения (при u (t ) В среде Matlab имеется символьная функция jacobian (.) , которая позволяет вычислять матрицу Якоби в символьном виде. Обычно, синтаксис команды вызова этой функции имеет вид: jacobian ( F , v) , где F - вектор – функция (в символьном виде), v - вектор переменных по F (i ) которым вычисляется матрица Якоби. Элемент (i, j ) результата определяется, как . v( j ) Пример. Вычислим матрицу Якоби векторной функции F ( x, y ) и x2  y2  2 . Файл E16_Jacob.   x(6  2 x  y) y(4  x  y)  в окрестности точек x1  y1  0 14 К сожалению, символьное вычисление матрицы Якоби в среде Matlab, не всегда возможно. В этом случае, для ее вычисления необходимо использовать приближенные методы, где вычисление производных заменяется их разностной аппроксимацией. Довольно часто, при линеаризации, возникают более сложные случаи. Предположим, что исследуемая система (объект) следует вдоль некой программной траектории xˆ (t ), uˆ (t ) , где xˆ (t )  F ( xˆ(t ), uˆ(t )) . Определим, как уже указывалось выше, x(t )  x(t )  xˆ(t ) , u(t )  u(t )  uˆ (t ) . Тогда мы можем записать уравнения: F x  xˆ (t ) F x  xˆ (t ) xˆ (t )  x(t )  F ( xˆ (t ), uˆ (t ))  |u uˆ (t ) x(t )  |u uˆ (t ) u . x u Так как xˆ (t )  F ( xˆ (t ), uˆ (t )) , то получим: F x  x F x  x x(t )  |u u x(t )  |u u u . x u F x  xˆ(t ) F x  xˆ (t ) |u uˆ(t ) nn и B(t )  |u uˆ (t ) nm зависят от В этом случае матрицы Якоби A(t )  x u параметра t (или от времени). То есть линеаризованная модель первого приближения не является LTI моделью. Такие модели называются линеаризованными моделями первого приближения относительно программной траектории xˆ (t ), uˆ (t ) . Основные модели представления детерминированных линейных систем. Представление линейных систем в виде моделей состояния. Представление детерминированной линейной системы в виде модели уравнений состояния можно записать, в общем случае, в виде: x  A(t ) x  B(t )u, y  C(t ) x  D(t )u , где x  n , u  m , y   p , A(t )  nn , B(t )  nm , C(t )   pn , D   pm 15 Модель состояния стационарной линейной системы имеет, соотвественно, форму x  Ax  Bu , y  Cx  Du , где x n , u   m , y   p , A nn , B   nm , C  pn , D  pm . Операторное представление линейных стационарных систем. Прямое преобразование Лапласа функции f (t ) определяется с помощью соотношения:  L{ f (t )}   e st f (t )dt  F ( s) , где обозначение L{ f (t )} используется для обозначения интеграла Лапласа для краткости. Функция: F (s) является функцией комплексной переменной s    i   . Функция f (t ) называется оригиналом, а функция F (s) - изображением. Преобразование Лапласа определено для функции f (t ) , удовлетворяющей следующим условиям: 1. Для значений 2. Для значений 3. Функция t 0 t 0 функция имеет значение функция f (t )  0 ; f (t ) имеет конечное число точек разрыва непрерывности; f (t ) ограничена в росте, при t   , экспонентой | f (t ) | M  et , где |  |  , M - неотрицательная константа Обратное преобразование Лапласа, при заданном изображении F (s) определяется следующим образом: где  i  1 f (t )  L {F ( s )}  F ( s )e st ds  2  i i  больше, чем реальная часть сингулярности изображения F (s) . 1 , Рассмотрим LTI (Linear Time Invariance) модель системы, которая описывается линейным дифференциальным уравнением n - го порядка с нулевыми начальными условиями: x ( n )  an 1 x ( n 1)  ...  a1 x  a0 x  bm u ( m )  bm 1u ( m 1)  ...  b1u   b0u ; x(0)  0 , x(0)  0 , …, x ( n 1) (0)  0 , u(0)  0 , u(0)  0 ,…, u ( m 1) (0)  0 . Применив преобразование Лапласа, получим, что: bm s m  bm 1s m1  ...  b0 X ( s)  n  U ( s) . s  an 1s n 1  ...  a0 Таким образом, каждой паре вход - выход {x(t ), u(t )} для LTI модели можно поставить в X ( s) соответствие пару {X (s),U (s)} на s – плоскости. При этом соотношение будет U ( s) определяться дробно – рациональной функцией от параметра s в виде: X ( s) bm s m  bm1s m1  ...  b0  n U ( s) s  an 1s n 1  ...  a0 которая называется передаточной функцией. Определение. , 16 Передаточной функцией системы управления называется отношение изображения X ( s ) выходного сигнала к изображению U ( s) входного сигнала при нулевом состоянии системы (нулевых начальных условиях). Для LTI модели системы передаточная функция представляет собой дробно – рациональную функцию параметра s в виде: bm s m  bm 1s m 1  ...  b0 W ( s)  n s  an 1s n 1  ...  a0 . Такое представление системы в среде Matlab называется tf – формой. Наивысший порядок n знаменателя передаточной функции определяет порядок системы. Для физически реализуемых систем порядок числителя m не должен превышать порядок знаменателя, то есть m  n . В этом случае передаточная функция называется правильной. Если выполняется соотношение m  n , то такая функция называется строго правильной. Величина называется относительной степенью передаточной функции. Модель системы в виде передаточной функции (tf – модель) реализуется в среде Matlab с помощью следующей команды: W  tf [num, den] , где num  [bm , bm1 ,..., b0 ] , den  [1, an1 , an2 ,..., a0 ] . Пример. Передаточную функцию (точнее передаточную матрицу) для систем с многими входами и выходами можно получить из уравнений линеаризованной модели первого приближения системы при нулевых начальных условиях. Пусть задано уравнения состояния системы вида: x  Ax  Bu , y  Cx  Du , где x n , u   m , y p . Такая форма записи системы в среде Matlab называется ss – моделью. Применяя преобразование Лапласа к левой и правой части уравнения найдем: Y (s)  (C  (s)  B  D)  U (s) ,  ( s )  ( sE  A) 1 . Отсюда получим, что передаточная матрица системы равна: W (s)  C  (s)  B  D . где: Анализ линейных стационарных систем при операторном представлении. Частотные характеристики Под термином «частотный выход системы» понимают установившееся значение выхода системы при подаче на ее вход синусоидального сигнала с различной частотой. Рассмотрим установившийся режим при синусоидальном входе для устойчивых, линейных, стационарных систем, которые описываются с помощью передаточной функции. Напомним, что устойчивые системы имеют все полюса, расположенные в левой части комплексной плоскости, то есть вещественные части полюсов являются отрицательными. Обозначим передаточную функции 17 W ( s ) , а, соответственно изображение входа через L{x(t )}  X (s) и выхода, как L{ y(t )}  Y (s) . Если вход x (t ) является синусоидальным сигналом некоторой частоты  , то системы, как выход системы в установившемся режиме будет также синусоидальным с той же частотой, но с другой амплитудой и фазой сдвига. Пусть входной сигнал имеет вид: x(t )  A  sin( t ) , где частота  измеряется в rad/sec. Представим передаточную функцию W (s) с простыми полюсами в виде: W ( s)  B( s ) B( s )  A( s) ( s  s1 )( s  s2 )  ...  ( s  sn ) Изображение выхода системы тогда равно: Y (s)  W (s)  X (s) . Декомпозицию изображения выхода, в виде разложения на простые дроби, можно представить как: Y ( s)  W ( s)  X ( s)  W ( s)  где c, c полюса d d  A c c    1  ...  n 2 s  s  i   s  i   s  s1 s  sn 2 , и di являются, в общем случае, комплексными константами. Для устойчивой системы s1 , s2 ,..., sn имеют отрицательные вещественные части. Поэтому, при больших значениях st s t s t , составляющие выхода e 1 , e 2 ,..., e n можно принять равными нулю. Таким образом, установившиеся значение можно принять равным: t yss (t )  c  e it  c  eit , где вычеты c и c могут быть вычислены следующим образом: 2 A A  W (i  ) , c  lim (s  i  )  W (s)  2   s  i  s  2 2i 2 A A  W (i  ) . c  lim ( s  i  )  W (s)  2  s i  s  2 2i Так как W (i  ) является комплексной величиной, можно, используя формулу Эйлера, записать: W (i  ) | W (i  ) | ei , где   arctan( Im(W (i  )) ). Re(W (i  )) С другой стороны имеет место равенство: W (i  ) | W (i  ) | e  i . Тогда установившийся выход равен: ei ( t )  ei ( t )  A | W (i  ) | sin( t  )  H sin( t  ) , 2i H  A | W (i  ) | . Отсюда можно сделать вывод, что установившийся выход yss (t )  A | W (i  ) | где переменная также представляет собой синусоидальный сигнал той же частоты, что и входной синусоидальный сигнал. Но амплитуда выходного сигнала и фаза (угол сдвига) будут, в общем случае, отличными от входного сигнала. На рисунке 2.3 показаны соотношения между входным сигналом и выходным сигналом системы. 18 Функция W (i  ) Рис. 2.3 называется частотной передаточной функцией или амплитудно-фазовой характеристикой (АФЧХ). Функция (АЧХ). Функция | W (i  ) | амплитудно-частотной характеристикой системы   ()  arg(W (i  ))  arctan( Im(W (i  )) ) Re(W (i  )) называется фазо-частотной характеристикой (ФЧХ) системы или частотным аргументом. Пример. Задана передаточная функция системы в форме Боде: W ( s)  k , где k , T Ts  1 вещественные константы. x(t )  A  sin( t ) . Тогда амплитудно-фазовая характеристика имеет k k вид: W (i  )  . Амплитудно-частотная характеристика равна: | W (i  ) | и 2 2 1 i T   1 T  фазо – частотная характеристика: ()   arctan(T  ) . Установившийся выходной сигнала будет k sin(  t  arctan(T  )) . описываться соотношением: yss (t )  A 1  T 2 2 Пусть задан синусоидальный вход: Так как передаточная функция W (s) определяется как:  W ( s)   e  st w(t )dt ,  где w(t ) - импульсная функция, то формально можно записать, что: W (i  )   eit w(t )dt . Таким образом, амплитудно-фазовая характеристика является Фурье преобразованием от импульсной функции. Так как имеет место равенство: e it  cos(  t )  i  sin(  t ) , то тогда можно записать следующие соотношения:   Re{W (i  )}   w(t ) cos(  t )dt , Im{W (i  )}    w(t )sin(  t )dt . Отсюда несложно получить: Re{W (i  )}  Re{W (i  )} , Im{W (i  )}   Im{W (i  )}. Таким образом, амплитудно-фазовая характеристика для частот из диапазона   0 является зеркальным отражением этой же характеристики для диапазона частот   0 относительно вещественной оси. 19 Если все полюса и нули передаточной функции W (s) имеют отрицательную вещественную часть, то такая передаточная функция называется минимально – фазовой. Если передаточная функция имеет хотя бы один нуль или один полюс с положительной вещественной частью, то такая система называется неминимально – фазовой. Если система является минимально – фазовой, тогда амплитудно – фазовая характеристика W (i  ) будет полностью определяться амплитудно – частотной характеристикой | W (i  ) | . Если система является неминимально – фазовой, то функция | W (i  ) | , так и W (i  ) будет определяться, как () . Правильная передаточная функция ( m  n ) будет ограничена по амплитуде с увеличением частоты, то есть | W ( i   ) |  . Строго правильная передаточная функция ( m  n ) будет стремиться к нулю при увеличении частоты, то есть | W (i  ) | 0 . Построение частотных характеристик. Амплитудно-частотная характеристика является комплексной функцией и зависит от параметра  . Имеется несколько вариантов графического представления такой характеристики, к которым относятся следующие: - Диаграммы Боде или логарифмические частотные характеристики; - Диаграмма Найквиста; - Диаграмма Николса. Диаграммы Боде. Логарифмические частотные характеристики. Диаграмма Боде состоит из двух графиков. Первый график представляет собой зависимость логарифма амплитудно – частотной характеристики (ЛАЧХ), а второй график отражает зависимость фазо-частотной характеристики (ЛФЧХ) от частоты входного синусоидального сигнала. Оба графика строятся в логарифмической шкале частот. Логарифмическая амплитудно-частотная характеристика строится в следующем масштабе: LogA()  20log10 | W (i  ) | . То есть единицей измерения функции Log ( A()) является децибел, которая имеет аббревиатуру dB. Таким образом, изменение амплитуды в 10 раз при изменении частоты в 10 раз будет приводить к изменению наклона логарифмической амплитудночастной характеристике на графике на 20dB / decade . Хотя график функции Log ( A()) невозможно построить во всем диапазоне частот, так как log10 0   , это не создает серьезных трудностей. Это обусловлено тем, что исследование реальных систем ведется на частотах отличных от нуля, то есть   0 . Передаточную функцию можно представить в виде произведения элементарных множителей вида: W ( s )  W j ( s ) , j Где W j (s) (Ts  1) , - передаточные функции некоторых элементарных звеньев, таких как: 1 1 2 , (T s  2T   1) , 2 Ts  1 T s  2T   1 и T1 s  1 T22 s  2T2   1 K , sv , 1 sv , . Если известны характеристики таких элементарных передаточных функций, то по ним легко определить и соответствующие характеристики передаточной функции W (s ) . Действительно, логарифмические амплитудно-частотную и фазо - частотную характеристики (ЛАЧХ и ЛФЧХ) можно вычислить с помощью следующих соотношений: 20 Log ( A())  20log10 | W j (i  ) | 20 log10 | W j (i  ) |   Log ( Ai ()) , j j j ()  arg(W (i  ))  arg(W j (i  ))   arg(W j (i  )) . j j Если ЛАЧХ Log ( A()) и ЛФЧХ () известны, то их можно соответствующим образом корректировать, добавляя в исходную передаточную функцию в качестве множителей передаточные функции некоторых элементарных звеньев. Построение диаграмм Боде в среде Matlab. Функция bode(.) средыMatlab позволяет строить диаграммы Боде для LTI моделей, заданных в виде передаточной функции. Данная функция имеет следующий синтаксис: bode(W ) - автоматическое построение диаграммы Боде; bode(W ,(m , M )) - построение диаграммы Боде в диапазоне частот (m , M ) ; bode(W , ) - построение точек диаграммы Боде для частот заданных вектором  ; [ A, , ]  bode(W ) - вывод данных диаграммы Боде в массив данных без построения графика; bode(W1 ,' ',W2 ,'*',W3 .' : r ') - построение единой диаграммы Боде для различных систем; bode( A, B, C, D) - построение диаграммы Боде для системы, заданной моделью состояния: y  Cx  D , где x n , y, u  Построение асимптотических логарифмических характеристик. Пусть передаточная функция системы представляет собой произведение передаточных функций элементарных звеньев, то есть: W ( s)  W j ( s) . Тогда ЛАЧХ и ЛФЧХ такой передаточной j функции можно представить в виде: Log ( A())  20  log 10 j | W j (i  ) |   Log ( Aj ()) , j ()   arg(W j (i  ))    j () . Для оценки правильности расчетов, а также для j j предварительного анализа такой системы, можно построить асимптотическую аппроксимацию ЛАЧХ и ЛФЧХ системы с передаточной функцией W ( s ) . Для этого можно использовать следующие правила: 1. Представить передаточную функцию W (s) в форме Боде; 2. Определить значение коэффициента усиления в логарифмической шкале также угловые (собственные) частоты элементарных звеньев упорядочить эти частоты по возрастающему порядку, то есть 3. Построить точку (  1, 20 log10 K ) 4. Провести через точку (  1, на диаграмме Боде 20 log10 K ) j  1 Tj 20log10 K ,а ; при этом 1  2  .... ; (log10 , dB) ; первую асимптотическую линию с наклоном 20v dB/decade, где v определяет порядок идеальных интегральных или дифференциальных звеньев sv до угловой частоты   1 . 21 5. Построить вторую асимптотическую линию для диапазона частот через точку (  1, 20 log10 K ) и имеющую наклон 20v 1    2 , проходящую плюс наклон асимптоты первого элементарного звена, имеющего угловую (собственную) частоту, совпадающую с 6. Построить асимптотические линии для последующих участков правилом пункта 5; наклон j   1 .  j 1     j , пользуясь - ой асимптоты определяется наклоном итоговой ЛАЧХ в точке  j 1 плюс наклон асимптоты j - го элементарного звена; 7. Построить последнюю асимптоту, проходящую через последнюю угловую (собственную) частоту, и продолжить ее до    . Диаграмма Найквиста. Диаграмма Найквиста для АФЧХ W (i  ) строится в координатах Re{W (i  )},Im{W (i)} при изменении частоты  от 0 до  ( в среде Matlab обычно принимается, что  изменяется от  до  ). График W (i  ) на диаграмме представляет кривую, которую описывает конечная точка вектора ( Re{W (i  )},Im{W (i)}) при изменении частоты . Таким образом, проекции точки на оси координат определяют вещественную и мнимую части АФЧХ W (i  ) . Построение диаграммы Найквиста в среде Matlab. Для построения диаграммы Найквиста в среде Matlab используется функция nyquist (.) . Команда ее выполнения может иметь следующий синтаксис: nyquist (W ) - автоматическое (с параметрами по умолчанию) построение диаграммы Найквиста; [Re,Im, ]  bode(W ) - вычисляет координат точек диаграммы Найквиста без построения графика; nyquist (num, den) - построение диаграммы Найквиста для передаточной функции numerator polynomial in s , заданной числителем и знаменателем; W ( s)  deno min ator polynomial in s nyquist (num, den, ) - построение диаграммы Найквиста для фиксированных частот, заданных вектором  ; 22 nyquist ( A, B, C, D, iu) - построение диаграммы Найквиста для модели состояния x  Ax  Bu ; y  Cx  D , где x n , y, iu  , параметр iu определяет вход системы. Анализ линейных систем при представлении в виде уравнений состояния. Рассмотрим автономную линейную систему следующего вида: x  Ax  Bu; y  Cx  Du; x n , u m , y  p , где A, B, C , D - действительные матрицы размерностей (n  n),(n  m),( p  n),( p  m) соответственно. Обычно предполагается, что начало координат x0  0 является положением равновесия для свободной, или неуправляемой системы, для которой u (t )  0 . При этом считается, что более общая система с ненулевым положением равновесия ( x , u ) во многих случаях может быть приведена к такому же виду с помощью параллельного переноса осей координат. Для такого приведения необходимо уметь решать линейное уравнение, которое определяет точку x  0 при значении u  0 , то есть уравнение вида Ax  b . Переход от уравнения состояния к передаточной функции. Пусть модель состояния системы представляется соотношениями: x  Ax  Bu; y  Cx  Du; x n , u m , y  p . Эквивалентная tf модель будет иметь матричную передаточную функцию вида: W ( s )  C  ( s  E  A) 1 B  D , где E - единичная диагональная матрица размера ( n  n ) . Для вычисления такой передаточной функции можно использовать метода Леверье – Фадеева. Обозначим характеристический полином системы в виде: P( s )  det( s  E  A)   n s n   n 1s n 1  ...   0 , где n  1 . Тогда обратная матрица (резольвента): ( s  E  A) 1 может быть записана в форме: Qn 1s n 1  Qn  2 s n  2  ...  Q0 ( s  E  A)  P( s ) 1 где Qn 1  E , . Метод Леверье – Фадеева позволяет одновременно вычислять коэффициенты характеристического полинома i , i  0,1, 2,..., n  1 и матрицы Q j , j  0,1, 2,..., n  2 используя следующие рекуррентные соотношения / 4 /: Qn 1  E ,  n  1 ; Fn 1  A  Qn 1 ,  n 1  tr ( Fn 1 ) / 1; Qn  2  Fn 1   n 1  E ; Fn  2  A  Qn  2 ,  n  2  tr ( Fn  2 ) / 2; Qn 3  Fn  2   n  2  E ......... F1  A  Q1 , 1  tr ( F1 ) / (n  1); Q0  F1  1  E ; F0  A  Q0 ,  0  tr ( F0 ) / n; F0   0  E  0 . 23 n Здесь: tr ( F )   Fii . Последнее соотношение может быть использовано для вычисления i 1 обратной матрицы A 1  Q0 0 . Преобразования математических моделей в среде Matlab. Matlab позволяет эффективно преобразовывать различные модели LTI систем . Команда [ A, B, C, D]  tf 2ss(num, den) позволяет преобразовать tf модель в модель состояния. При этом, передаточная функция задается в виде векторов коэффициентов числителя и знаменателя: W ( s)  numerator polynomial in s . deno min ator polynomial in s Обратное преобразование от модели состояния (ss модели) к модели передаточной функции (tf модели) осуществляется командой [num, den]  ss2tf ( A, B, C, D, iu) , где параметр iu должен указываться для систем, имеющих более, чем один вход. Если система имеет только один вход, то команда может иметь синтаксис: [num, den]  ss2tf ( A, B, C, D) . Внутренняя устойчивость. Устойчивость по Ляпунову. Понятие внутренней устойчивости или устойчивости по Ляпунову основывается на оценке количественного поведения систем при нулевом входном сигнале, то есть исследуется устойчивость собственного движения системы исключительно от ненулевых начальных условий. Рассмотрим нелинейную систему: x  F (t , x) , F (t , 0)  0 , x(t0 )  x 0 , x n , которая имеет изолированную точку равновесия в начале координат. Будем считать, что функция F (.) является непрерывной и имеет непрерывные частные производные t  t0 . Тогда решение системы x  x(t )  n F , k  1, 2,..., n xk можно представить как траекторию в при n - мерном пространстве в зависимости от независимого параметра t . Определение. Решение   (t ) системы x  F (t , x) называется устойчивым по Ляпунову при t   (или просто устойчивым), если: 1. Для любого   0 существует такая величина 2. (, t0 )  0 , что все решения x (t ) системы x  F (t , x) (включая   (t ) ) удовлетворяют условию || x(t0 )  (t0 ) ||  ; Все решения определены на всем интервале [t0 , ) и удовлетворяют неравенству || x(t )  (t ) ||  для t  [t0 , ) . Интуитивно понятно, что если решение   (t ) является устойчивым, то другие  x (t ) , которые были достаточно близки к решению  в момент t0 , будут для t  [t0 , ) лежать в некоторой заданной узкой трубке, построенной вокруг решения   (t ) решения x 24 F (t , 0)  0 то тривиальное решение (t )  0 (называемое также невозмущенным движением системы) будет устойчивым, если для любого t0 существует такая величина (, t0 )  0 , из неравенства || x(t0 ) ||  будет вытекать следующее соотношение || x(t ) ||  для t  [t0 , ) . Следует понимать, однако, что если решение   (t ) является устойчивым, то это не значит, что все решения будут ограничены. Таким образом, если все решения системы x  F (t , x) В частности, когда являются ограниченными, то это не значит, что они устойчивы, и наоборот. Если можно выбрать величину   0 такую, что она не будет зависеть от параметра t 0 , то есть   () , тогда такое решение   (t ) называется равномерно устойчивым. Определение. Решение   (t ) называется неустойчивым по Ляпунову, если существует такое решение x  x (t ) и величины   0 ,   0 такие, что || x (t0 )  (t0 ) ||  , и найдется такой момент t1  t1 ()  t0 , при котором будет выполняться неравенство || x (t1 )  (t1 ) ||  . В частности, когда F (t ,0)  0 , то тривиальное решение  (t )  0 будет неустойчивым, если существует такое решение x  x (t ) и величины   0 ,   0 такие, что || x(t 0 ) ||  , и найдется такой момент t1  t1 ( )  t 0 , при котором будет выполняться неравенство Определение. Решение   (t ) называется асимптотически устойчивым при || x(t1 ) ||  . t   , если выполняются следующие условия: - Решение (t ) является устойчивым; -   (t0 )  0 , t0 , что все решения x  x (t ) , которые удовлетворяют условию || x(t0 )  (t0 ) ||  , имеют свойство lim || x(t )  (t ) || 0 Существует такая величина t  25 В частности, если решение условия || x(t0 ) ||  (t )  0 является асимптотически устойчивым, то при выполнении будет выполняться свойство lim || x(t ) || 0 . t  Устойчивость LTI моделей систем. Рассмотрим однородную LTI модель системы: x  Ax , x (t 0 )  x 0 , x   n . Общее решение системы можно записать в виде x(t )  e A(t  t 0 ) 0 x . Построим преобразование системы в эквивалентную жорданову каноническую форму с помощью преобразования z  Tx , где T  R n n несингулярная матрица. Тогда можно записать следующее соотношение: A  T 1 ( J1 (1 )  ...  J s (s ))T T , где 1 ,...,  s , s  n - собственные значения матрицы A с учетом их кратности. Обозначим через m1 ,..., ms порядки жордановых клеток (ячеек). Отсюда можно записать следующее соотношение: exp{( t  t0 ) I p ( p )}  e  p (t  t 0 ) [ Em p  (t  t0 ) (t  t0 ) (t  t0 ) I1, p  I 2, p  ...  I m 1, p ] , 1! 2! (m p  1)! p где I k , p  ( I p ) k , I p - единичная косая матрица размера m p  m p , E m p - единичная диагональная матрица размера m p  m p . Теорема. Если однородная LTI модель системы является x  Ax устойчивой по Ляпунову, то тогда все собственные значения 1 ,..., n матрицы A  R n n имеют неотрицательные вещественные части, то есть Re(k ))  0; k  1,2,..., n . Кроме того, алгебраическая кратность всех собственных значений, которые имеют нулевую вещественную часть, должна совпадать с их геометрической кратностью. Отсюда вытекает, что если LTI модель устойчива, то она будет равномерно устойчива. Действительно, если система x  Ax устойчива, то имеет место оценка || e At || M , где M является константой t  t0 . Тогда можно записать неравенство: || x |||| e A(t  t 0 ) ||  || x 0 || M || x 0 || . Если    , тогда получим: || x ||  . То есть   0 и не зависит от t 0 . Таким образом выбрать: || x || M тривиальное решение x(t )  0 является равномерно устойчивым. Теорема. Однородная LTI модель системы x  Ax будет асимптотически глобально устойчивой тогда, когда все собственные значения матрицы A  R n n имеют отрицательные вещественные части, то есть Re(k )  0; k  1,2,..., n . Очевидно, что все решения LTI моделей систем x  Ax являются или устойчивыми, или неустойчивыми одновременно. К сожалению, такого нельзя сказать о нелинейных системах, где некоторые решения могут быть устойчивыми, а другие неустойчивыми. Устойчивость по Ляпунову неоднородных LTI моделей систем. Рассмотрим неоднородную LTI модель системы: x  Ax  Bu , x  R n , u  Rm , 26 где вход u  u(t ) является ограниченной функцией, и который воздействует на систему на u (t ), t  [t0 , T ] , || u(t ) || um   . Общее  0, t  T ограниченном интервале [t 0 , T ] , T   , то есть u (t )   решение системы имеет вид: x  e t x   e A(t  ) u ( )d . Предположим, что соответствующая A( t  t 0 ) 0 t0 однородная система x  Ax является асимптотически устойчивой, то есть   max(Re (k ))  0 , где k k , k  1,2,..., n собственные значения матрицы A . Тогда можно записать следующие соотношения: || x |||| e t A( t  t 0 ) 0 x ||  ||  e A(t  ) u ( )d || e  (t  t 0 )  || x || u m  || e A(t  ) || d . Так как t0 t T A(t  ) || d |  e (t  ) d | et | [  || e t0 t0 T t0 e t 0  e T  ] | , и   0 , то отсюда получим || x(t ) || 0 при t   . Можно также показать, что если собственные значения с максимальными вещественными частями имеют одинаковую алгебраическую и геометрическую кратности, то если однородная система является устойчивой, то и неоднородная соответствующая система с ограниченным входом будет также устойчивой. Таким образом, устойчивость неоднородной LTI модели системы с ограниченным входным воздействием будет полностью определяться устойчивостью соответствующей однородной LTI модели. Оценка устойчивости автономной нелинейной системы по устойчивости линейной системы первого приближения. Пусть задана нелинейная автономная система x  F (x) , F (0)  0 . Построим в окрестности начала координат LTI модель первого приближения вида x  Ax и найдем связь между ее устойчивостью и устойчивостью исходной нелинейной системы. Теорема (1-я теорема Ляпунова об устойчивости). Если все собственные значения матрицы A линеаризованной стационарной системы первого приближения лежат в левой части комплексной плоскости, то невозмущенное движение соответствующей (то есть тривиальное решение) соответствующей (порождающей) нелинейной автономной системы x  F (x) будет устойчивым. Если, по меньшей мере, имеется одно собственное значение матрицы A лежит в правой части комплексной плоскости, то тогда невозмущенное движение соответствующей нелинейной автономной системы x  F (x) будет неустойчивым. В случае, когда имеются собственные значения, лежащие на мнимой оси, но нет, ни одного собственного значения, лежащего в правой полуплоскости, называется критическим. В критическом случае нельзя оценить устойчивость исходной нелинейной системы на основе математической модели линеаризованной стационарной системы первого приближения и требуется проведение дополнительного исследования. Внешняя устойчивость. BIBO устойчивость. Внешняя устойчивость еще называется устойчивостью по входу – выходу и определяет зависимость выхода системы, при воздействии на нее ограниченного входа. То есть, это принципиально отличается от внутренней устойчивости или устойчивости по Ляпунову, где в 27 качестве возмущения рассматривается только отклонение в начальных условиях. Так как при внешней устойчивости исследуется ограниченность выхода при подаче на вход ограниченного сигнала (Bounded -Input- Bounded-Output или BIBO), то часто такой тип устойчивости, особенно в англоязычной литературе называют BIBO устойчивостью. Однако, хотя внутренняя и внешняя устойчивости по основным признакам отличаются, они являются в определенном смысле связанными, что будет показано позднее. Говорят, что вход u (t ) является ограниченным, если u (t ) не стремится по амплитуде к   или, что эквивалентно, существует такая константа u m , что || u(t ) || um   для t  0 . Будем говорить, что система называется BIBO устойчивой, если для каждого ограниченного входа существует свой ограниченный выход. Данный тип устойчивости определяется при нулевом состоянии системы, то есть, когда все начальные состояния ее координат равны нулю. Другими словами, система будет BIBO устойчива, если существует такая положительная константа  , что при подаче на вход системы ограниченного сигнала u (t ) , при нулевом состоянии системы, выход системы будет удовлетворять неравенству sup || y (t ) ||   sup || u (t ) || . t 0 t 0 BIBO устойчивость SISO LTI систем. Рассмотрим SISO LTI систему, заданную передаточной функцией W (s ) . Тогда изображение выхода системы будет равно Y ( s)  W (s)U ( s) , где U (s)  L(u(t )) , | u(t ) | um   , или во временной t  области y (t )  w(t   )u ( )d . Теорема. t  SISO LTI система y (t )  w(t   )u ( )d будет BIBO устойчивой, тогда и только тогда, когда ее импульсная функция w(t ) будет абсолютно интегрируема на интервале [0, ) , или   | w(t ) | dt  M   , где M - некоторая положительная константа. Рассмотрим свойства BIBO устойчивости для SISO LTI систем. Теорема. Если система с импульсной функцией w(t ) является BIBO устойчивой, то тогда при t   , получим: 1. Если вход равен u (t )  a для t  0 то выход определяется соотношением W (0)  a ; 2. Если вход равен u(t )  sin( 0t ) для t  0 , то выход определяется соотношением | W (i  0 ) |  sin( 0t  arg W (i  0 )) , где: W ( s )    w(t )e  st dt . Теорема. SISO система, которая задается правильной дробно – рациональной передаточной функцией W (s) будет BIBO устойчивой, тогда и только тогда, когда каждый полюс W (s) имеет отрицательную 28 вещественную часть или, что эквивалентно, лежит в левой части комплексной плоскости (s – плоскости). BIBO устойчивость MIMO LTI систем. Теорема. MIMO LTI система, с правильной дробно - рациональной передаточной матрицей G( s)  ( g ij ( s)) , i  1,2,..., p; j  1,2,..., m , является BIBO устойчивой, тогда и только тогда, когда каждый полюс каждого элемента gij (s) имеет отрицательную вещественную часть. Рассмотрим связь между полюсами передаточной матрицы G(s) и собственными значениями C  [ Adj( sE  A)]B  D , то очевидно, что каждый полюс G(s) является det( sE  A) собственным значением матрицы A . То есть, если каждое собственное значение имеет матрицы A . Так как G ( s)  отрицательную вещественную часть, то и каждый полюс будет иметь отрицательную вещественную часть и, следовательно, система будет BIBO устойчивой. Однако, при вычислении передаточной матрицы возможно сокращение полюсов и нулей, в том числе и положительных. Поэтому, если исходная матрица A имеет собственные значения с нулевой или даже положительной частями, то система x  Ax  Bu , y  Cx  Du может быть BIBO устойчивой . То есть, используя понятие BIBO устойчивости, ничего нельзя сказать о внутренней устойчивости при не нулевых начальных условиях и нулевом входе. Управляемость линейных систем. Интуитивно, понятие управляемости обозначает возможность перевода состояния системы из одной точки пространства состояний в другую с помощью допустимого управления. Событие (t0 , x0 ) линейной системы x  A(t ) x  B(t )u называется управляемым относительно события (t1 , x1 ) , если существует момент времени t1  t 0 и вход u (t ) , который определен на конечном интервале [t 0 , t1 ] и, обеспечивает перемещение системы из события (t0 , x0 ) в событие (t1 , x1 ) . Для линейных систем управляемость по состоянию, обычно, осуществляется по отношению к точке x1  0 , так как если систему управляема по отношению к точке начала координат, то она будет управляема и относительно любой точки x1 . Теорема. Линейная система x  A(t ) x  B(t )u является управляемой, если существует управление u (t ) , которое позволяет переместить систему из точки x(t0 )  x0 в точку x(t1 )  x1 , когда вектор x0  (t0 , t1 ) x1 принадлежащую области линейного преобразования: t1 W (t0 , t1 )   (t0 , t ) B(t ) BT (t )T (t0 , t )dt , где матричная функция  (t , t0 )   1 (t0 , t ) iявляется t0  переходной функцией однородной линейной системы: x  A(t ) . Кроме того, если x является  T T  решением уравнения: W (t0 , t1 ) x  x0   (t0 , t1 ) x1 , то тогда вход u (t )   B (t ) (t0 , t1 ) x является возможным управлением, которое позволяет переместить систему из точки x(t0 )  x0 в точку x(t1 )  x1 за интервал времени [t 0 , t1 ] . 29 Следствие. Если в некоторые моменты времени t 0 и t1 , матрица W (t0 , t1) имеет ранг равный n , то тогда система x  A(t ) x  B(t )u является управляемой. Таким образом, управляемая система может быть преобразована в начало координат из любой другой точки за любой малый интервал времени. Это будет возможно, если пара матриц ( A, B) полностью известна и удовлетворяет условию W (t0 , t1 )  0 . Кроме того, нет других ограничений на сигнал управления. Однако, в этом случае, сигнал управления может быть очень большим. Рассмотрим теперь случай LTI модели системы x  Ax  Bu , x  R n , u  R m . Если LTI модель системы является управляемой, то тогда говорят, что пара матриц ( A, B) является управляемой. Для простоты, в дальнейшем примем, что t 0  0 . Теорема. Следующие утверждения являются эквивалентными: 1. Матричная пара ( A, B) является управляемой; t t T T 2. Матрица W (t )  e A BB T e A  d  e A(t  ) BB T e A (t  ) d является несингулярной для   любого t  0 ; 3. Матрица управляемости F  ( B, AB,..., An1B) , размера (n  nm) имеет ранг n (то есть полный построчный ранг); 4. Если все собственные значения матрицы A имеют отрицательные вещественные части, T T то тогда единственное решение WG уравнения: AWG  WG A   BB является положительно определенным. Решение WG называется Грамианом управляемости и  T задается соотношением: WG  e A BB T e A  d .  Следствие. Пара матриц ( A, B), A  R n n , B  R n m будет управляемой тогда и только тогда, когда матрица Fn  r  ( B, AB ,..., A n  r B ) имеет ранг равный n , где r  rank (B) . Пример. Рассмотрим систему стабилизации платформы / 2 /. 30   0.5 0   0.5   , B    . Начальное состояние  1  0  1  Система описывается уравнением is x  Ax  Bu , где A   платформы x1 (0)  10 , x2 (0)  1 . Найти управление, которое приведет платформу в состояние равновесия за 2 секунды. Ранг матрицы управляемости равен: rankF  0.5  0.25    2 , то есть система  rank ( B, AB)  rank   1   1 управляема. Вычислим Грамиан: 2   0.5  0.5   e  0.5 0       . 5 1    0 e  t  1    0.2162 0.3167  0    . Тогда управление имеет вид: d   e  t   0.3167 0.4908   e 0.5( 2  t ) 0  1  e 1 1  10  At1 T AT (t1  t ) 1    . u (t )   B e W (t1 )(e x0  x1 )  0.5 1 W (2)  2  e  ( 2  t )    0 e   1 e W (t  2)    0 0 Рассмотрим, изменяется ли управляемость LTI модели, если провести линейное x  Tx , где T является несингулярной матрицей. преобразование координат системы ~ Теорема. Управляемость LTI моделей является инвариантом по отношению ко всем невырожденным эквивалентным преобразованиям. Лемма. SISO LTI модель, которая задается дифференциальным уравнением вида: x ( n )  qn 1 x ( n 1)  ...  q0 x  p0u , p0  0 является управляемой. Область пространства состояний, где все точки являются управляемыми, называется областью управляемости LTI модели системы. Если линейная система является управляемой, то это значит, что она полностью управляема во всем пространстве состояний. Если линейная система является неуправляемой, то есть ранг матрицы управляемости rankF  r , 0  r  n , то тогда говорят, что система является частично управляемой. Для частично управляемой системы существует подпространство (область) управляемости, которое определяется базисом из независимых столбцов матрицы F . Функции оценки управляемости LTI моделей среды Matlab. Функция ctrb(.) позволяет строить матрицу управляемости LTI модели. Синтаксис команды ctrb(sys) позволяет формировать матрицу управляемости F на основе дескриптора (указателя) системы sys . Синтаксис ctrb( A, B) формирует матрицу F на основе ссылки на пару матриц ( A, B) . Функция gram(.) вычисляет Грамиан управляемости путем решения матричного уравнения Ляпунова. Наблюдаемость линейных систем. Понятие наблюдаемости можно охарактеризовать, как возможность оценки состояния системы по данным наблюдений за выходом. Проблема наблюдаемости системы отличается от проблемы ее идентификации. Например, для LTI моделей, проблема идентификации сводится к нахождению оценок матриц A, B, C, D на основе данных измерений входа и выхода системы. Задача 31 же наблюдения сводится к оценке текущего (или начального) состояния системы по измерениям входа и выхода. LTI модель состояния: x  Ax  Bu , y  Cx  Du называется наблюдаемой, если для любого неизвестного начального состояния x (t 0 ) , существует такой конечный момент времени t1  t 0 , что по результатам измерений входа u (t ) и выхода y(t ) на интервале [t 0 , t1 ] можно получить оценку состояния x (t 0 ) . В противном случае, модель называется ненаблюдаемой. Выход модели: x  Ax  Bu , y  Cx  Du , который обусловлен влиянием начального ненулевого состояния x (t 0 ) и входным воздействием u (t ) , имеет вид: t y (t )  Ce A(t  t 0 ) x(0)  C  e A(t  ) Bu ( )d  Du (t ) . Для простоты, будем в дальнейшем полагать, что: t0 t 0  0 . При этом предполагается, что вход u (t ) и выход y(t ) известны, а начальное состояние x(0) y (t )  Ce At x(0) , которую также можно определить с неизвестно. Введем обозначение функции ~ t помощью выражения: ~ y (t )  y (t )  C  e A(t  ) Bu ( )d  Du (t ) . Теорема. Следующие утверждения эквивалентны: 1. Пара матриц ( A, C ) является наблюдаемой, где A  R n n , C  R p n ; t T 2. Матрица Wo (t )  e A t C T Ce At dt является несингулярной для любого t  0 ;   C     CA  3. Матрица наблюдаемости: H   , размера ( pn  n) , имеет ранг равный n (полный ...   n 1   CA    ранг по столбцам). 4. Если все собственные значения матрицы A имеют отрицательные вещественные части, то матричное уравнение Ляпунова AT Wo  Wo A  C T C имеет единственное решение Wo , которое является положительно определенным. Решение Wo называется  T Грамианом управляемости и имеет вид: Wo  e A t C T Ce At dt .  Следствие. Пара матриц ( A, C ) является наблюдаемой, тогда и только тогда, когда T T T T nr T T C ) имеет ранг равный n , где r  rank (C ) . матрица H n  r  (C , A C ,..., ( A ) Рассмотрим линейное эквивалентное преобразование ~ x  Qx модели состояния, где Q несингулярная матрица. Теорема. Свойство наблюдаемости является инвариантным по отношению к невырожденному линейному преобразованию. 32 Функции оценки наблюдаемость LTI систем в среде Matlab. Функция obsv(.) позволяет вычислить матрицу наблюдаемости LTI модели системы. Синтаксис команды obsv(sys) формирует матрицу наблюдаемости H на основе дескриптора (указателя) sys системы. Команда obsv( A, C ) позволяет вычислить матрицу H по паре матриц ( A, C ) системы. Функция gram(.) позволяет вычислить Грамиан наблюдаемости, путем решения соответствующего матричного уравнения Ляпунова. Синтаксис команды имеет вид: gram(sys, ' o' ) . Принцип дуальности (двойственности). Построим следующую модель состояния вида: xˆ  AT xˆ  C T uˆ, yˆ  BT xˆ  DT uˆ , где: xˆ  R n , uˆ  R p , yˆ  R m . Такая модель системы называется двойственной по отношению к системе: x  Ax  Bu , y  Cx  Du , x  R n , u  R m , y  R p . Теорема (принцип двойственности). Пара матриц ( A, B) является управляемой, тогда и только тогда, когда пара матриц ( AT , B T ) является наблюдаемой и, соответственно, пара матриц ( A, C ) будет наблюдаемой, тогда и только тогда, когда пара ( AT , C T ) является управляемой. Декомпозиция линейных стационарных систем по Калману Рассмотрим LTI модель: x  Ax  Bu , y  Cx  Du , x  R n , u  R m , y  R p . Пусть задано x  Px , где матрица P является несингулярной. Тогда, преобразованная система преобразование ~ ~ ~ ~ ~ ~ ~ ~ ~ x  A~ x  Bu, y  C~ x  Du , где A  PAP1 , B  PB, C  CP 1 , D  D . имеет вид: ~ Теорема (декомпозиция управляемости). Пусть задана LTI модель: x  Ax  Bu , y  Cx  Du , x  R n , u  R m , y  R p где матрица управляемости имеет ранг: rank ( F )  rank ( B, AB,..., An 1B)  r  n . Тогда существует такая 1 матрица Q  P  (q1 , q2 ,..., qr , qr 1 ,..., qn ) , первые r столбцов которой состоят из r независимых столбцов матрицы F , а оставшиеся (n  r ) столбцов формируются произвольно, x  Px таким образом, чтобы матрица P была несингулярной. При этом преобразование ~ позволяет записать преобразованную модель состояния системы в виде:  Ac  xc   xc   x   A  xuc   Bu   0  uc   A12   xc   Bc  y  (C , C )   xc   Du ,    u, x  c uc Auc   xuc   0   uc  ~ ~ ~ rr ( n  r ) ( n  r ) где матрица управляемости F  PF , Ac  R , Auc  R , а подмодель состояния: ~ ~ ~ ~ xc  Ac ~ xc  Bc u , ~ y  Cc ~ xc  Du , ~ xc  R r , u  R m , y  R p является управляемой и имеет передаточную функцию, совпадающую с передаточной функцией исходной модели состояния системы. Пример. 33 0 1 0 0 1 1 2 0 1    0       Задана модель состояния x  Ax  Bu , y  Cx  Du , где A   1 0 0  , B   0  , C T   0  , D  0 .    0 0 0   Вычислим матрицу управляемости F  (b, Ab, A b)   0 0 0  . Получим rank ( F )  1 . Построим 1 0 0   2 0 1 1   матрицу P   0 1 0  , det( P)  1 , где первый столбец P является первым столбцом матрицы F , а другие 1 0 0   столбцы выбираются из условия, что det( P)  0 . Тогда преобразованная система может быть записана в ~ ~ x   A ~ ~ x  A~ x  B u или  c    c форме: ~ ~    xuc   0  ~ ~ 1 xc   Bc  0 2 2 ~ A12  ~ ~   1         u , где , A  PAP  0 0 1 B  PB   0  , ~  ~      Auc  xuc   0   0 0 1 0   ~ ~ ~ ~ C  C  P 1  0 0 1 , Ac  0 , Bc  1 , Auc   0 1  . 1 0 Теорема (декомпозиция наблюдаемости). Пусть задана модель состояния: x  Ax  Bu , y  Cx  Du , x  R n , u  R m , y  R p , матрица  C     CA   l  n . Тогда существует такая наблюдаемости которой имеет ранг: rank ( H )  rank  ...   n 1   CA     p1   ...    матрица P   pl  , где первые l строк состоят из l независимых транспонированных  ...   pn  столбцов матрицы H , а оставшиеся строки матрицы P выбираются таким образом, чтобы x  Px позволяет записать она была несингулярной. Тогда линейное преобразование ~ эквивалентную систему в форме: ~ ~ ~  Ao ~ x  ~ xo   Bo  x  ~ ~ 0  ~ ~  xo     o   A ~ c   B     u~ ~  ~    ~ u , y  (Co ,0) ~   Du , где ~    xuc   xuo   xuo   A21 Au 0  xuo   Buo  ~ ~ Ao  R l l , Auo  R ( n  l ) ( n  l ) , и подмодель состояния: ~ ~ ~ ~ xo  Ao ~ xo  Bo u , ~ y  Co ~ xo  Du , ~ xc  R l , u  R m , y  R p является наблюдаемой и имеет передаточную функцию, равную передаточной функции исходной системы. Комбинируя декомпозицию управляемости и декомпозицию наблюдаемости можно перейти к общей декомпозиции моделей состояния по Калману. Теорема (Декомпозиция LTI моделей по Калману). Каждая модель состояния может быть преобразована некоторым линейным преобразованием к следующей форме: 34 ~  ~ xc 0   Aco   ~ ~ xcuo   A21  ~   xuco   0 ~     xucuo   0 ~ A13 ~ A23 ~ Auco ~ A43 ~ xc 0   Bco  0  ~  ~  ~  xcuo   Bcuo  A24  ~ ~ ~ x  Du , где вектор ~  u , y  Cco 0 Cuco 0  ~ xco  ~   0  xuco   0  ~ xucuo   0  Aucuo  ~ x является управляемым и наблюдаемым; вектор ~ является управляемым, но не наблюдаемым; ~ Acuo   cuo xuco является наблюдаемым, но не управляемым, и вектор ~ xucuo является не вектор ~ управляемым и не наблюдаемым. Более того, исходная модель состояния эквивалентна, при нулевых начальных условиях, управляемой и наблюдаемой ее редуцированной части: ~ ~ ~ ~ xco  Aco ~ xco  Bco u , ~ y  Cco ~ xco  Du , и имеет передаточную функцию равную: ~ ~ ~ ~ W ( s )  Cco ( sE  Aco ) 1 Bco  D . Результаты данной теоремы могут быть проиллюстрированы в графическом виде. Из рисунка видно, что передаточная функция описывает только часть существующей модели состояния системы. Функции декомпозиции LTI моделей в Matlab. Функция ctrbf (.) позволяет осуществлять преобразование моделей состояния, которые задаются тройкой матриц ( A, B, C ) , в форму декомпозиции управляемости. Синтаксис команды [ Ac, Bc, Cc,T , k )  ctrbf ( A, B, C ) позволяет построить декомпозицию модели, заданную ~ A 0 0 ~ ~ матрицами Ac, Bc, Cc , представляемых в виде: Ac   ~uc ~  , Bc   ~  , Cc  Cuc Cc . A   Bc   21 Ac  Матрица преобразования T определяется как ортогональная несингулярная матрица. Вектор k показывает, какое число итераций было проведено алгоритмов декомпозиции на каждом   шаге. Величина sum(k ) указывает на число независимых столбцов в матрице управляемости, а также определяет размерность управляемой части модели состояния. Функция obsvf (.) позволяет преобразовать модель состояния, заданную тройкой матриц ( A, B, C ) , в форму декомпозиции наблюдаемости. Синтаксис команды [ Ao, Bo, Co, T , k )  ctrbf ( A, B, C ) позволяет вычислять матрицы Ao, Bo, Co , формы ~ ~ ~  Buo  A12  ~   Bo  , ~  ~  , Co  0 Co . Матрица   0 Ao   Bo   преобразования T ищется как ортогональная и несингулярная. Вектор k определяет число итераций A декомпозиции наблюдаемости в виде: Ao   uo   35 алгоритма на каждом шаге. Величина sum(k ) определяет число независимых строк в матрице наблюдаемости, а также размерность наблюдаемой части преобразованной модели. Минимальная реализация SISO систем. Понятие грубость линейных систем. Нам уже известно, что система, заданная передаточной функцией, может иметь несколько моделей состояния, заданных разными уравнениями одной размерности. Одновременно, передаточная функция является в ряде случаев избыточным описанием, так имеет сокращаемые нули и полюса. Это порождает проблему поиска такой реализации передаточной функции, чтобы получить минимальную размерность уравнения состояния среди всех возможных. В данном разделе мы рассмотрим задачу поиска минимальной реализации для SISO LTI систем, которая тесно связана с проблемой грубости линейных систем. Пусть система с одним входом и одним выходом задана правильной передаточной функцией вида: P(s) pm s m  pm 1s m 1  ...  p0  ,m  n. Q( s ) s n  qn 1s n 1  ...  q0 Будем называть передаточную функцию W (s ) несократимой, нет нулей и полюсов, которые можно W ( s)  взаимно сократить, или, что эквивалентно, числитель и знаменатель не имеют общих корней или множителей. n Рассмотрим теперь модель состояния SISO системы: x  Ax  Bu , y  Cx , где x  R , u, y  R . Эта система имеет передаточную функцию W ( s )  C ( sE n  n  A) 1 B . Предположим теперь, что один или более нулей и полюсов можно взаимно сократить, то есть сократить общие множители в числителе ~ P (s) ~ ~ ~  deg ree( P P(s ) и знаменателе Q(s) . Тогда можно записать: W ( s )  ~ , где m (s))  m , Q(s) ~ n~  deg ree(Q(s))  n . ~ Ясно, что модель состояния, построенная по передаточной функции W ( s ) , имеет ~ , в то время, как исходная система имеет порядок n  n~ , то есть не является размерность n ~ ~ ~ минимальной. Обозначим через ( A, B, C ) матрицы модели состояния, построенной по ~ ~ 1 ~ ~ ~ передаточной функции W ( s ) : W ( s )  C ( sE n~  n~  A) B . Возникает вопрос, является ли такая реализация системы минимальной? Теорема. Пусть задана передаточная функция W ( s)  pm s m  pm 1s m 1  ...  p0 где n  deg ree(Q(s)) s n  qn 1s n 1  ...  q0 степень знаменателя передаточной функции. Тогда n - мерная модель состояния, реализованная по передаточной функции W (s ) , будет минимальной тогда и только тогда, когда передаточная функция W (s) является несократимой. Пример. Задана модель состояния системы x  Ax  Bu ,  0 1 0   1  ;   6  11  6    y  Cx , где x  R 3 , u, y  R , A   0 1  0  6( s  1)   ; T   . Отсюда найдем передаточную функцию . Очевидно, что W ( s)  B   6  C  0 ( s  1 )( s  2 )( s  3 ) 0   30      36 несократимая передаточная функция равна ~ ~ W (s)  ~ системы будет модель состояния: x  A x  B u , 6 . Тогда минимальной реализацией ( s  2)( s  3) 1  ~ 0 ~  0 ~  , B    , y  Cx , где x  R 2 , u, y  R , A     6  5 6 1 C T    .  0 Пример.   4  1.25 1.5   2     0  , B  0 , Пусть задана система x  Ax  Bu , y  Cx , где x  R , u, y  R , A   4  0 0 1 0      0    6( s  1) T C   0.75  . Отсюда найдем передаточную функцию W ( s)  , то есть ( s  1)( s  2)( s  3)   0.75    3 ~ W (s)  ~ 6 ~ ~ 2 . Тогда минимальная реализация имеет вид: is: x  A x  B u , y  Cx , where x  R , ( s  2)( s  3) 1  ~  0 T 1 ~  0  , B    , C    . u, y  R , A     6  5 6  0 Таким образом, из анализа примеров, мы видим, что две различные системы имеют одну и ту же минимальную реализацию. Рассмотрим теперь поведение исходных систем и их минимальной реализации. Все траектории системы первого примера и ее минимальной реализации стремятся к точке равновесия, которая для обеих систем имеет один и тот же тип: устойчивый узел. Следовательно, минимальная реализация имеет такое же поведение, что и исходная система. Во втором примере часть траекторий исходной системы стремится к точке равновесия, а другая часть траекторий удаляется от точки равновесия. То есть точка равновесия имеет тип: седло. У системы минимальной реализации все траектории стремятся к точке равновесия типа устойчивый узел. То есть поведение исходной системы и ее минимальной реализации различно. Следовательно, если для дальнейшего проектирования выбрать минимальную реализацию, то при проектировании можно совершить ошибку. Малейшее различие в параметрах исходной системы приведет к невозможности сокращения нулей и полюсов, лежащих в правой части комплексной плоскости. Это приведет к тому, что система станет неустойчивой. Отсюда следует, что при построении несократимой передаточной функции можно взаимно сокращать только нули и полюса, лежащие в левой части комплексной плоскости. Понятие грубости линейных систем. Интуитивно понятно, что грубость систем, как математических моделей, заключается в том, что небольшие изменения их параметров должны приводить к небольшим изменениям в их поведении. Впервые понятие грубости систем было введено А. Андроновым и Л. Понтрягиным. Физическое обоснование грубости систем заключается в следующем: если математическая модель известна приблизительно, а такое имеет место в большинстве практических случаев, то дискуссия о влиянии на поведение системы малых изменений ее параметров является бессмысленной. Математический аспект проблемы грубости заключается в том, что малые изменения параметров системы не должны изменять топологический портрет ее траекторий, то измененная система и 37 исходная должны быть эквивалентными, за исключением особых критических случаев, требующих дополнительного исследования (так называемые точки бифуркации параметров). Основной вывод из приведенных рассуждений является следующим: проектируемая система должна быть грубой, то есть малые изменения параметров не должны изменять ее качественное поведение и переводить в другой топологический тип, за исключением критических случаев, связанных с потерей устойчивости. Управляемость, наблюдаемость и минимальная реализация MIMO систем. Рассмотрим SISO LTI модель системы, которая задается строго правильной передаточной функцией: W ( s)  D( s ) , deg ree( D( s))  deg ree( A( s)) . Предположим, что нули k , k  1,2,..., n A( s) полинома A(s ) являются простыми и вещественными. То есть можно записать соотношение: W ( s)  n D (k ) D ( k ) . Очевидно, что величины bk  являются константами. A(k ) k )( s  k )  A( k 1 Построим по передаточной функции модель состояния, используя преобразование Лурье: n xk  k xk  bk u, k  1,2,..., n , y   xk . Тогда модель состояния примет вид: x  Ax  Bu , y  Cx , где k 1  b1   1 ... 0      B   ...  , A   0 ... 0  , C  1 ... 1 . Вычислим матрицу управляемости такой модели: b   0 ...   n  n   b1  b F  ( B, AB..., A n 1 B )   2  . b  n 1 1  n 1  2 det( F )   bk  det k 1 . . 1  n  b11 b2 2 . bn n ... b11n 1  ... b2 n2 1   . Отсюда найдем: ... .  ... bn nn 1  ... 1n 1  ... n2 1  . ... .  ... nn 1  Второй множитель, в выражении определителя матрицы управляемости, является определителем Вандермонда, который не равен нулю. Тогда определитель матрицы F будет равен нулю, когда существует хотя бы один коэффициент bk  0 . НО это будет возможно, когда нуль числителя D(s) равен s  k . То есть, это возможно, когда передаточная функция имеет сокращаемые нуль и полюс. В этом случае, координата xk модели состояния будет не зависеть от входа u . Таким образом, передаточная функция, имеющая сокращаемые нули и полюса, не может описывать полностью модель состояния. Рассмотрим LTI модель x  Ax  Bu , y  Cx  Du , x  R n , u  R m , y  R p . Примем, для простоты, D  0 . Пусть передаточная функция системы W ( s)  C ( sE  A) 1 B является несокращаемой, то есть полиномы числителя и знаменателя не содержат общих множителей / 5, 6/. Тогда можно сформулировать следующие условия полной наблюдаемости и управляемости. 38 Теорема. Модель состояния x  Ax  Bu , y  Cx, x  R n , u  R m , y  R p является полностью управляемой и наблюдаемой, тогда и только тогда, когда передаточная функция W ( s)  C ( sE  A) 1 B является несокращаемой. Если начальное состояние модели состояния является нулевым, то выход y(t ) системы ( A, B, C, D) может быть получен, на основании декомпозиции по Калману, в виде: t ~ ~ ~ y (t )  Cco  e Aco (t  ) Bco u ( )d  Du (t ) . Данное выражение совпадает с выходом полностью t0 управляемой и наблюдаемой части системы. Для простоты примем, что D  0 . Тогда ~ ~ ~ подсистема ( Aco , Bco , Cco ) будет совпадать с минимальной реализацией исходной системы. Минимальная реализация модели всегда является полностью управляемой и наблюдаемой, и имеет несокращаемую передаточную функцию. Алгоритм получения минимальной реализации модели состояния может быть следующим (здесь формы декомпозиции управляемости и наблюдаемости приняты совпадающими с соответствующими формами среды Matlab): 1 1. Найти матрицу преобразования Tc , позволяющую выделить управляемую и ~  Auc ~ 1 неуправляемую части (подсистемы): A  Tc ATc   ~ A  21 ~ ~ ~ C  CTc  Cuc Cc .   0 0 ~ 1 ~  , B  Tc B   ~  , Ac   Bc  ~ ~ ~ 2. Найти матрицу преобразования To такую, что в управляемой подсистеме ( Ac , Bc , Cc ) ~ ~  Acuo Ac,12  ~ 1 ~ позволяет выделить наблюдаемую часть (подсистему): Ao  To AcTo   ~ ,  0 Aco   ~ ~ ~ ~ ~ B  ~ Bo  To1Bc   ~cuo  , Co  CcTo  0 Cco .  Bco  0   En  rank ( A~ ) c  . Тогда матрица T 1  Tˆo1Tc1 будет определять 3. Построить матрицу: Tˆo1   1   To   ~ ~ ~ преобразование исходной системы в минимальную ее реализацию в виде ( Aco , Bco , Cco ) .   С помощью линейного преобразования с матрицей T исходная система преобразуется в форму ~  Auc ~ декомпозиции по Калману: z   A21   0 ~ Acuo 0   0  ~  ~  ~ Ac,12  z   Bcuo u , y  Cuc ~  ~ Aco   Bco   Пример. Задана модель состояния x  Ax  Bu , y  Cx  Du , x  R , u  R, y  R , где: 4  ~ 0 Cco z  Du . 39 5   4 A   0  0  4     7 0 0   2 , B , C  2 2  2 1 , D  0 .  0 0 4 2      1  0  2 6    8 Выделить управляемую и наблюдаемую подсистему (минимальную реализацию) . Ранг матрицы управляемости F  ( B, AB, A2 B, A3 B) равен 3. Ранг матрицы наблюдаемости H  (C T , AT C T , ( AT ) 2 C T , ( AT )3 C T )T равен 3. Передаточная функция исходной системы в zpk 10( s  2.6)( s  3)( s  4) представлении имеет вид: W ( s )  . Таким образом, порядок минимальной ( s  1)( s  2)( s  3)( s  4) 10( s  2.6) реализации системы равен 2. Передаточная функция минимальной реализации равна Wmr ( s )  . ( s  1)( s  2) В среде Matlab имеется функция minreal(.) , которая позволяет формировать минимальную реализацию LTI моделей систем. С помощью этой функции из системы удаляются неуправляемая и ненаблюдаемая части. Анализ линейных нестационарных систем. Асимптотическая устойчивость линейной стационарной системы полностью определяется корнями характеристического многочлена ее коэффициентов. Поэтому, естественно, попытаться связать асимптотическую устойчивость неавтономной системы x  A(t ) x с корнями  j (t ) уравнения: det[ A(t )  E ]  0 . В первую очередь желательно ответить на вопрос: гарантирует ли неравенство: max Re( j (t ))   , j где   0 - положительная постоянная для любого t [t0 , ) , асимптотическую устойчивость системы? Ответ на этот вопрос отрицательный. Пример.   x  1  2cos 4t 2  2sin 4t  x . Собственными значениями ее 2  2sin 4t 1  2cos 4t характеристического многочлена являются функции 1 (t )  2 (t )  1 . Однако эта система обладает Пусть задана система: неограниченным решением 1 (t )  et sin 2t ,  2 (t )  et cos 2t . Таким образом, в общем случае, неравенство max Re( j (t ))   для корней j характеристического многочлена не гарантирует асимптотическую устойчивость неавтономной системы. Переходная матрица автономной линейной системы вида x  Ax , имеет следующий вид: X (t )  exp( At ) . При этом, элементы фундаментальной матрицы, нормированной в точке t  0 , представляют собой линейные комбинации функций p j (t ) exp(  j t ) , где p j (t ) - некоторые 40 многочлены,  j - собственные числа матрицы A . В терминах этих чисел (собственных значений), как известно, можно полностью характеризовать свойство асимптотической устойчивости системы  x  Ax . Оказывается, как показал А.М. Ляпунов, со всякой неавтономной линейной дифференциальной системой можно связать конечный набор действительных чисел, которые, как и в автономном случае, позволяют исследовать ее асимптотическую устойчивость. Пусть f (t ) - числовая функция, заданная и непрерывная на полупрямой R . Число ( или символы   и  ) определяемое формулой: ____ 1 t  t [ f ]  lim ln | f (t ) | , называется характеристическим показателем или показателем Ляпунова функции f (t ) . Здесь под ____ обозначением lim  (t ) понимается верхний предел функции t   (t ) при t   . Рассмотрим теперь однородную неавтономную систему x  A(t ) x с непрерывной на матрицей A(t )  R n*n R . Теорема Ляпунова (о характеристических показателях решений линейной неавтономной системы) . Если матрица A(t ) ограничена, то есть || A(t ) || c  const   , то каждое ненулевое решение системы x  A(t ) x имеет конечный характеристический показатель. Определение. Множество всех характеристических показателей ненулевых решений системы x  A(t ) x называется ее спектром. Теорема. Спектр линейной системы x  A(t ) x с ограниченной матрицей A(t ) состоит из не более, чем n различных элементов. Используя понятие спектра, можно достаточно детально описать структуру общего решения уравнения x  A(t ) x . Пусть совокупность функций x1 (t ), x2 (t ),..., xn (t ) образует фундаментальную систему, а числа  1 ,  2 ,...,  n -полный спектр рассматриваемого уравнения, то есть  [ x j ]   j , j  1,2,..., n . Положим  j  e  j t x j (t ) . Отсюда получим  [ j ]   j   [ x j (t )]  0 . Отсюда, общее решение уравнения x  A(t ) x может быть представлено в виде n  c j j (t )e  jt , где j 1 c j - произвольные постоянные, а функции  j имеют нулевые характеристические показатели. Отсюда сразу же следует, что если полный спектр системы отрицательный, то есть  j  0 , то такая система асимптотически устойчива. Теорема. Неравенство Важевского. Для любого решения системы x  A(t ) x , где A(t )  C , при t 0 где t t t0 t0  t   справедливо неравенство: || x(t 0 ) || exp{   ( )d } || x(t ) |||| x(t 0 ) || exp{  ( )d } ,  (t ) и  (t ) наименьший и наибольший характеристический корни эрмитово- симметризованной матрицы: A (t )  H 1 [ A(t )  A* (t )] . 2 41 Следствие 1. Для асимптотической устойчивости системы x  A(t ) x достаточно выполнения условия  (t )  h  0 , где h - положительная постоянная. Нетрудно убедиться, что спектр системы x  A(t ) x отрицательный, если квадратичная T H форма x A (t ) x определенно отрицательна, то есть справедливо условие x T A H (t ) x   | x | 2 , где   const  0 , а A H (t ) в случае действительной матрицы A(t ) имеет вид A (t )  [ A(t )  A (t )] . Данное условие имеет место, когда главные миноры  k (t ) матрицы H T A H (t ) отделены от нуля, то есть det{ k (t )}    const  0 и чередуют знаки. Причем первый из определителей меньше нуля. Критерий неустойчивости неавтономных систем. Пусть задана система x  A(t ) x . Тогда 1t lim  tr ( A( ))d  0 t  t t0 ____ из условия следует неустойчивость ненулевых решений системы. Методы и алгоритмы синтеза линейных стационарных систем Частотные методы синтеза систем. Можно сформулировать следующую схему решения задачи синтеза в частотном диапазоне. 1. Построить ЛАЧХ Lnc () и ЛФЧХ  nc ( ) (диаграмму Боде) некомпенсированной разомкнутой системы. 2. На основе требований к точности и качеству замкнутой системы построить ЛАЧХ Ld ( ) и ЛФЧХ  d ( ) желаемой компенсированной разомкнутой системы. Получит, методом графического вычитания, частотные характеристики компенсатора: Lmc ()  Ld ()  Lnc () , mc ()  d ()  nc () . 3. 4. Построить, на основе анализа частотных характеристик компенсатора: Lmc (), mc () , передаточную функцию компенсатора. 5. Вычислить передаточную функцию замкнутой системы при включении в ее состав спроектированного компенсатора. Оценить устойчивость, точность и качество замкнутой системы. 6. Сравнить полученные оценки с заданными требованиями к замкнутой системы. Если требования выполняются, то приступить к проектированию физически реализуемого компенсатора. Если требования не выполняются, то вернутся к шагу 2 алгоритма и откорректировать частотные характеристик желаемой разомкнутой системы. Как видно, из приведенных выше рассуждений и примеров, основной проблемой синтеза систем в частотном диапазоне является построение желаемых частотных характеристик компенсированной разомкнутой системы, которая будет удовлетворять заданным требованиям по точности и качеству. Укажем несколько практических рекомендаций, которые облегчают процесс построения желаемых частотных характеристик. Прежде всего, необходимо определить частотный диапазон работы замкнутой системы путем анализа частотных характеристик входа, выхода и шумов. Затем выделить низко – частотный, средне – частотный и высоко - частотный диапазоны. Низко – частотный диапазон определяется минимальной угловой частотой l наиболее медленного звена, входящего в состав некомпенсированной разомкнутой системы. ЛАЧХ разомкнутой системы имеет наклон 0 dB или кратный  20 v dB/decade, где параметр v определяет 42 порядок астатизма системы. Амплитуда ЛАЧХ на частоте   1 определяет желаемый коэффициент усиления системы. Средне – частотный диапазон определяется частотами в интервале c1    c2 , который содержит частоту среза   c ; c [c1, c2 ] желаемой ЛАЧХ разомкнутой системы. В соответствии с критерием Найквиста частота среза должна удовлетворять соотношению c  b , где частота  b определяется условием  d (b )  180 0 . Обычно предполагается, что в окрестности частоты:   c желаемые частотные характеристики должны совпадать с соответствующими частотными характеристиками устойчивых звеньев первого или второго порядка, то есть наклон желаемой ЛАЧХ должен составлять -20 dB/decade или -40dB/decade. Величина  c выбирается из условия обеспечения требуемого времени t s затухания переходного процесса и максимального P перерегулирования  . Для этого можно воспользоваться таблицей значений функции t s c ( max ) , P0 где Pmax , P0 - максимальное и начальное значения амплитудно – частотной характеристики Pc ( ) замкнутой системы. P Используя заданное значение перерегулирования  r , можно по таблице найти величину max , а P0 затем найти произведение t sc . Зная требуемое значение t s , отсюда легко найти частоту среза  c желаемой разомкнутой системы. Затем, по оценкам запасов устойчивости по амплитуде и фазе, можно определить коэффициенты усиления: Ld (c1 ), Ld (c 2 ) . Частоты  c1 ,  c 2 выбираются таким образом, чтобы c1  10  Ld ( c1 ) 20  c и c 2  10  Ld ( c 2 ) 20  c , а также выполнялись неравенства:  c 2  10 , 2  c  4 .  c1  c1 Участок стыковки ЛАЧХ низко – частотного и средне – частотного диапазонов (то есть, когда l    c1 ) выполняется в виде прямых линий с наклоном -40 dB/decade или -60 dB/decade. Высоко – частотный диапазон желаемой ЛАЧХ (то есть при   c 2 ) выполняется с наклоном превышающим -40 dB/decade. Методы компенсации и желаемой передаточной функции. Метод желаемого размещения полюсов (root – locus method) является одним из основных способов улучшения динамики замкнутой системы во временной области, позволяющий достигать 43 требуемых прямых показателей качества. При этом используются различные схемы включения компенсаторов (контроллеров) в систему. На рисунках a), b) показаны последовательные схемы включения компенсаторов. Параллельная схема включения компенсатора приведена на рисунке c). Пример. Рассмотрим структурную схему, приведенную на рисунке a). Пусть передаточная функция объекта 2 управления равна W p ( s)  2 . Спроектируем контроллер, который будет гарантировать s  3s  1 время затухания переходного процесса не более t s  3 секунд, перерегулирование   0% и ошибку в установившемся режиме e()  0 . k s 2  as  b Для этого будем использовать PID контроллер (регулятор): Wcm ( s )  k p  k d s  i  k d , s где a  kp kd ,b  s ki . Выберем в качестве желаемой передаточной функции, удовлетворяющей kd 1 . Замкнутая передаточная функция системы с PID s 1 Wcm ( s)W p ( s) указанным выше требованиям: Wd ( s )  контроллером равна: Wc ( s)  1  Wcm ( s)W p ( s) . Выберем в качестве компенсатора динамики объекта контроллер с передаточной функцией вида: s 2  as  b  s 2  3s  1, i.e. a  передаточная функция замкнутой системы примет вид: Wc ( s)  kp kd 2k d  1 s  2k d  3, b  1 2k d s 1 ki  1 . Тогда kd  Wd ( s) . Отсюда найдем, что: kd  0.5 . Следовательно, получим: k p  1.5, ki  1 и, окончательно, получим, что передаточная функция контроллера равна: Wcm ( s )  1.5  0.5s  1 . s 44 Рассмотрим случай, когда желаемая передаточная функция последовательного контроллера Wd (s) замкнутой системы известна и удовлетворяет заданным требованиям. Тогда передаточную функцию компенсатора (контроллера) Wcm (s ) можно найти из уравнения: Wcm ( s)  Wcm ( s)W p ( s) 1  Wcm ( s)W p ( s)  Wd ( s) . Решение Wcm (s ) данного уравнения имеет вид: Wd ( s) m( s ) 1 . Перепишем передаточную функцию объекта в виде: W p ( s)  . Тогда W p ( s) (1  Wd ( s)) p( s) можно записать соотношение: Wcm ( s )  p( s ) Wd ( s ) . Чтобы такая функция была физически m( s ) (1  Wd ( s )) w(t ), t  0 , то есть:  0, t  0 реализована, она должна удовлетворять соотношениям: L1{Wcm ( s )}   Wcm ( s)  rl s l  ...  r0 s k  qk 1s k 1  ...  q0 , l  k . При этом, замкнутая система должна быть грубой, то есть при компенсации не происходило сокращение положительных полюсов и нулей на s – плоскости. Детализируем передаточную функцию объекта в виде: W p ( s)  m( s) m  ( s)m  ( s) , где p  ( s), m  ( s)    p( s) p ( s) p ( s ) полиномы, имеющие нули с отрицательной вещественной частью, а p  ( s), m  ( s) - полиномы, имеющие нули с положительной вещественной частью. Такое представление передаточной функции называется ее факторизацией / 2 /. Тогда передаточную функцию компенсатора можно записать в виде: Wcm ( s)  p  ( s) p  ( s) Wd (s) . Отсюда, чтобы нули и полюса с положительной m (s)m (s) (1  Wd (s))   вещественной частью взаимно не сокращались, желаемая передаточная функция должна удовлетворять условиям: Wd ( s )  m  ( s) K (s) p  ( s ) D( s ) , 1  Wd ( s)  , где полином N (s) - является N (s) N ( s) знаменателем передаточной функции Wd (s) , а K ( s), D( s) - некоторые полиномы, которые должны удовлетворять уравнению: m  ( s) K ( s)  p  ( s) D( s)  N ( s) . Данное уравнение будет разрешимо относительно K ( s), D( s) , когда степень полинома N (s) будет меньше, чем deg ree( K (s))  deg ree( D(s))  2 , то есть: deg ree( N (s))  deg ree( K (s))  deg ree( D(s))  1 (a). Если передаточная функция объекта имеет, в своем составе, v p интегральных звеньев, а передаточная функция замкнутой системы должна иметь астатизм v - го порядка (тип системы будет равен v ), то передаточная функция компенсатора будет содержать множитель s 1  Wd ( s)  p  ( s ) D( s ) s N ( s) примет вид: v v p , или: Wcm ( s )  m (s) K (s)  p  (s) D(s)s p  (s) K (s) m ( s) D( s) s v  v p v v p  v v p , то есть: . Тогда полиномиальное уравнение  N (s) . Тогда условие физической реализуемости компенсатора можно записать, как: deg ree( p  ( s))  deg ree( K ( s))  deg ree(m ( s))  deg ree( D( s))  v  v p (b). 45 Так как степень числителя передаточной функции Wd (s) не может превышать степень знаменателя, то тогда функция (1  Wd (s)) должна иметь числитель и знаменатель одного порядка. То есть получим следующее условие: deg ree( N ( s))  deg ree( p  ( s))  deg ree( D( s))  v  v p (c ). Отсюда, из анализа условий (a ), (b ), (c ) получим неравенство / 2,3 /: deg ree( N ( s))  deg ree( p( s))  deg ree( p  ( s))  deg ree(m ( s))  v  v p  1 . Таким образом, алгоритм определения передаточной функции компенсатора можно записать в следующем виде: m  ( s) K (s) 1. Определить желаемую передаточную функцию Wd ( s )  замкнутой системы со N ( s) знаменателем N (s) . m  ( s )m  ( s ) 2. Провести факторизацию передаточной функции объекта: W p ( s)  3. Записать условия (a ), (b ), (c ) и выбрать соответствующие минимальные степени p  ( s) p  ( s) . nK , nD полиномов K ( s), D( s) . 4. Записать выражения полиномов K ( s), D( s) с неопределенными коэффициентами и построить полиномиальное уравнение: 5. Решить уравнение m  m ( s) K (s)  p  (s) D(s)  N (s) . ( s ) K ( s )  p  ( s ) D ( s )  N ( s ) путем нахождения соответствующих коэффициентов K ( s), D( s) . 6. Задать передаточную функцию компенсатора в виде выражения: p  ( s) K ( s) Wcm ( s)   m ( s ) D( s ) s v  v p . Пример. Пусть передаточная функция объекта управления равна: W p ( s )  1 . Провести синтез ( s  1)( s  1) компенсатора, который будет обеспечивать замкнутой системе астатизм первого порядка и характеристический полином (знаменатель) равный: N ( s)  ( s  1) 3 . Таким образом, получим, что: v  1 . Факторизованная передаточная функция объекта имеет: m  ( s)  m  ( s)  1 , p  ( s)  s  1, p  ( s)  s  1 . Условия физической реализации последовательного компенсатора имеют вид: deg ree( N ( s))  3 , deg ree(m  ( s))  deg ree(m  ( s))  0 , deg ree( p  ( s))  deg ree( p  ( s))  1 . Тогда условия (a ), (b), (c ) можно записать в виде: (a ) 3  nK  nD  1 , (b ) 1  nK  nD  1 , (c ) 3  1  nD  1 . Отсюда найдем: nD  deg ree( D( s))  1; nK  deg ree( K ( s))  1 . Следовательно: D(s)  d1s  d0 ; K (s)  k1s  k0 . 3 2 Полиномиальное уравнение тогда имеет вид: k1s  k 0  ( s  1)( d1s  d 0 ) s  s  3s  3s  1 . Отсюда найдем: d1  1, d0  4 , k1  7, k0  1 . Следовательно, передаточная функция компенсатора может 46 быть записана в форме: Wcm ( s )  p  (s)  K (s) m ( s) D( s) s v v p  ( s  1)(7 s  1) , а желаемую передаточную ( s  4) s m ( s ) K ( s ) 7 s  1 функцию можно записать, как: Wd ( s)   N ( s) ( s  1)3 . Как мы уже обсуждали, имеется тесная связь между прямыми показателями качества замкнутой системы и расположением полюсов и нулей ее передаточной функции, которые определяются видом желаемой передаточной функции. Ограничения на прямые показатели качества обычно характеризуют не один, какой – то, конкретный переходный процесс, а целый класс временных ответов системы на стандартный вход / 4 /. Возникает вопрос, как в этом случае определить такой класс желаемых передаточных функций. Теорема. Пусть SISO LTI модель системы задана нормализованной передаточной функцией вида: (bm s m  ...  b1s  b0 ) ( p1 p2 ... pn ) Wn ( s )  , m  n , где полюса pi являются вещественными и n b0  (s  pi ) i 1 представлены в упорядоченном виде: 0  p1  p2  ...  pn . Обозначим переходную функцию такой системы Wn (s ) , как xnm (t , p1 , p2 ,..., p n ) . Соответственно, обозначим переходную функцию системы Wdn ( s)  (bm s m  ...  b1s  b0 ) ( p n ) m (t , p ) . Тогда существует такое вещественной , как x dn n b0 ( s  p) m (t , p(t )) для любого t  [0, ] , число p  0 , что будет выполняться равенство: xnm (t , p1 , p2 ,..., pn )  xdn при условии, что: p1  p  p n . Существует и более общий результат. Пусть задана SISO LTI модель системы, заданная передаточной функцией: Wn ( s)  bm s m  ...  b1s  b0 a0 , m  n , которая имеет s n  an 1s n 1  ...  a0 b0 полюса: pi , i  1,2,..., n с отрицательной вещественной частью, такими, что: Re( pi ) [ pmin , pmax ] , где: pmin  min Re( pi )) , pmax  max Re( pi )) , и переходной функцией xnm (t , p1 , p2 ,..., pn ) . Обозначим i m i через: x dn (t , p ) переходную функцию системы с передаточной функцией с кратными полюсами: 47  bm s m  ...  b0 | p |2 q , n  2q   ( s  p ) q ( s  p ) q b0 Wdn ( s )   , такими, что: pr [ pmin , pmax ] , bm s m  ...  b0 pr  | p |2 q  , n  2q  1  ( s  p )( s  p) q ( s  p ) q b0 r  Re( p) [ pmin , pmax ] , Im( p) [ f max , f max ] , где f max  max Im( pi )) . Тогда существуют вещественное i число pr  0 , и комплексное число p, Re( p)  0 такие, что будет выполняться равенство: m xnm (t , p1 , p2 ,..., pn )  xdn (t , pr , p) для любого t  [0, ] , и pr [ pmin , pmax ] , Re( p) [ pmin , pmax ] , Im( p) [ f max , f max ] . Если определить связь между величинами: Re( p) , и Im( p) как: определяет сектор с корневым углом:   d , d   2 Re( p)  tan( d ) , где угол  d Im( p) , то тогда можно задать только два параметра: Re( p) [ pmin , pmax ] , pr [ pmin , pmax ] , для случая, когда: n  2q  1 , или один параметр вида: Re( p) [ pmin , pmax ] , когда n  2q . Здесь угол  d определяет колебательность системы. Синтез систем методом обратной связи по состоянию. Пусть задана SISO LTI система: x  Ax  Bu ; y  Cx; x  R n ; u  R, y  R , и управление имеет вид: u(t )  g  Kx , где g  R задающее воздействие, а вектор K имеет размер (1  n) . Тогда замкнутая система будет описываться уравнением: x  ( A  BK ) x  Bg ; y  Cx / 3 /. Теорема. Пусть SISO LTI модель системы x  Ax  Bu ; y  Cx; x  R n ; u  R, y  R является полностью n n 1  ...   0 с вещественными управляемой и, задан желаемый полином вида  d ( )     n 1 коэффициентами. Тогда существует вектор K  R1 n обратной связи по состоянию такой, что модель состояния x  ( A  BK ) x; y  Cx замкнутой системы имеет характеристический полином, совпадающий с желаемым полиномом, то есть det(E  A  BK )  d ( ) . Лемма. Пусть SISO LTI модель системы: x  Ax  Bu ; y  Cx; x  R n ; u  R, y  R является полностью n n 1  ...   0 с вещественными управляемой и, задан желаемый полином вида:  d ( )     n 1 коэффициентами. Тогда существует вектор K  R1 n обратной связи по состоянию такой, что 1 передаточная функция замкнутой системы Wc ( s )  C ( sE  A  B K ) B имеет знаменатель, 48 совпадающий с желаемым полиномом  d (s ) , а числитель совпадает с числителем передаточной функции Wo ( s )  C ( sE  A ) 1 B разомкнутой системы. Таким образом, использование обратной связи по состоянию позволяет управлять расположением полюсов, но не нулей, передаточной функции замкнутой системы. Алгоритм построения обратной связи по состоянию можно сформулировать в следующем виде. 1. Проверить управляемость разомкнутой системы путем вычисления ранга матрицы управляемости: Fo  ( B, AB ,..., A n 1 B) . 2. Если разомкнутая система управляема, то определить характеристический полином  A ( )  n  an 1n 1  ...  a0 , матрицы A, B канонической формы Фробениуса - Калмана, матрицу Fo  ( B , A B ,..., A n 1B ) и передаточную функцию Wo ( s)  C ( sE  A) 1 B  3. 4. p s n 1  ...  p0 P( s ) .  n n 1 n 1  A (s) s  an 1s  ...  a0 1 Найти матрицу преобразования T  Fo Fo . Записать желаемый характеристический полином замкнутой системы:  d ( )  n   n 1n 1  ...   0 . 5. Определить элементы вектора обратной связи K для канонической формы, на основе соотношений: k j   j  a j ; i  1,2,.., n . 6. Вычислить вектор обратной связи для исходной системы: K  KT . Записать передаточную функцию замкнутой системы в форме: Wc ( s )  C ( sE  A  B K ) 1 B  pn s n 1  pn 1s n  2  ...  p1 .  d ( s) Как мы видим, определение вектора K обратной связи сводится к поиску решения уравнения: det(E  A  BK )  d () . Лемма (формула Аккермана). Если пара матриц ( A, B) является полностью управляемой, то решение уравнения det(sE  A  BK )  d (s) определяется следующим соотношением: K  [0,0,..., 0,1][ B, AB ,..., A n 1B]1 d ( A) . MIMO LTI модель замкнутой системы с обратной связью по состоянию. Рассмотрим MIMO LTI модель системы: x  Ax  Bu ; y  Cx; x  R n ; u  R m , y  R p , m  1, p  1 . Теорема. Если пара матриц ( A, B) является полностью управляемой и, задан желаемый полином:  d ( s )  s n   n 1s n 1  ...   0 с вещественными коэффициентами, то тогда существует обратная связь по состоянию: u(t )   Kx(t ) , где K  R m n - вещественная постоянная матрица, которая, для замкнутой системы: x  ( A  BK ) x  Bg ; y  Cx , позволяет получить характеристический полином равный желаемому, то есть: det(sE  ( A  BK ))  d (s) . Функции среды Matlab для синтеза обратной связи по состоянию. 49 В среде Matlab имеются две функции ac ker(.) и place(.) , которые позволяют рассчитывать параметры обратной связи по состоянию. Команда K  ac ker( A, B, p) предназначена для вычисления вектора обратной связи K дляSISO LTI полностью управляемых моделей, задаваемых парой матриц ( A, B) . Желаемый характеристический полином задается множеством желаемых собственных значений матрицы ( A  BK ) , которое задается вектором - строкой p . Функция place(...) предназначена для вычисления матрицы K обратной связи по состоянию для MIMO LTI полностью управляемых моделей систем. Команда K  place( A, B, p) , где пара матриц ( A, B) является полностью управляемой, а вектор – строка p , определяет желаемые собственные значения матрицы ( A  BK ) , позволяет вычислять матрицу K . Следует отметить, что вектор p не может содержать кратные собственные значения (функция, в этом случае, формирует сообщение об ошибке). LTI модели систем, стабилизируемые обратной связью по состоянию. Если модель состояния является полностью управляемой, то с помощью обратной связи по состоянию можно расположить все собственные значения матрицы замкнутой системы в заданной области. Теперь рассмотрим случай, когда модель состояния не является полностью управляемой. Рассмотрим систему: x  Ax  Bu ; y  Cx; x  R n ; u  R m , y  R p , которая имеет матрицу управляемости F  ( B, AB,..., An 1B), rank ( F )  r  n . Тогда, используя декомпозицию по Калману,  zc   Ac     zuc   0 можно уравнение состояния записать в эквивалентной форме:  A12  zc   Bc      u , Auc  zuc   0  r (n  r ) ; Ac  R r  r ; A12  R r  ( n  r ) ; Auc  R ( n  r ) ( n  r ) ; Bc  R ( r  m) . Так как пара ( Ac , Bc ) где z c  R ; zuc  R является управляемой, то существует такой закон управления u(t )  Kc zc , что матрица ( Ac  Bc Kc ) замкнутой управляемой части системы, будет иметь желаемое распределение собственных значений. Тем не менее, необходимо, чтобы вся замкнутая система обладала внутренней устойчивостью. То есть, матрица Auc должна быть гурвицевой, или характеристический полином det( sE  Auc ) должен иметь нули, имеющие отрицательные вещественные части. Определение. MIMO LTI модель системы: x  Ax  Bu ; y  Cx; x  R n ; u  R m , y  R p , называется стабилизируемой, если существует закон управления u(t )   Kx(t ) такой, что система является асимптотические устойчивой. Теорема. Если MIMO LTI модель системы не является полностью управляемой, то она будет стабилизируемой, тогда и только тогда, когда матрица Auc неуправляемой части системы будет гурвицевой, то ее собственные значения будут иметь отрицательную вещественную часть. Напомним, что матрица обратной связи для исходной модели системы определяется с помощью соотношения: K  KP. 50 Синтез LTI SISO систем с обратной связью. Рассмотрим, следующую LTI SISO модель состояния разомкнутой системы вида: x  Ax  Bu ; y  Cx; x  R n , y, u  R . Как известно, если пара ( A, B) является полностью управляемой, то существует такой вектор K  R1 n обратной связи по состоянию, что замкнутая система: x  ( A  BK ) x  Bg , где g  R - задающее воздействие, будет иметь желаемый характеристический полином  d (s ) . Кроме того, если пара ( A, C ) будет полностью наблюдаемой, то передаточная функция замкнутой системы является несокращаемой и имеет вид: Wc ( s )  C ( sE  ( A  BK )) 1 B  Wo ( s)  C ( sE  A) 1 B  P( s) , где полином P(s ) равен числителю передаточной функции  d (s) P( s ) разомкнутой системы. det( sE  A) Таким образом, чтобы переходной процесс в замкнутой системе удовлетворял заданным требованиям качества, необходимо исследовать передаточную функцию Wc ( s )  P(s) , где  d (s) полином P(s ) уже известен. Поэтому, для выбора желаемой области полюсов можно использовать метод желаемой передаточной функции с кратными вещественными полюсами, которую можно записать в форме: Wd ( s)  P( s ) n ( s   ) n P(0) . Отсюда, нетрудно определить возможный интервал   [ min,  max] . Если выбрать некоторую внутреннюю точку интервала *  [ min,  max] , то тогда желаемый полином будет определяться, как  d ( s )  ( s   ) n . Данное условие возможно расширить, если выбрать более общий вид желаемой передаточной с кратными полюсами. После нахождения желаемого полинома  d (s ) можно найти вектор коэффициентов обратной связи путем использования функций ac ker(.) или place(.) среды Matlab. Синтез матрицы усилений замкнутой системы. Рассмотрим обратную связь по состоянию для LTI MIMO системы в форме: u   Kx  Hg , где K матрица обратной связи соответствующей размерности, H - матрица усилений размера (m p) и имеющая ранг rank ( H )  p , и g  R p вектор входных задающих воздействий. LTI модель разомкнутой системы имеет вид: x  Ax  Bu ; y  Cx; x  R n , u  R m , y  R p . Тогда уравнение замкнутой системы можно записать в виде: x  ( A  BK ) x  BHg ; y  Cx . Передаточная функция 1 замкнутой системы равна: Wc ( s )  C ( sE  ( A  BK )) BH . Заметим, что обеспечить требования к установившемуся значению каждого выхода можно только тогда, когда размерность задающего вектора входных воздействий g имеет размерность равную p / 2 /. Поэтому, обычно, должно выполняться ограничение m  p , то есть в разомкнутой системе число управляющих входов должно быть больше, чем число измеряемых выходов. Будем также предполагать, что LTI модель разомкнутой системы будет полностью управляемой и наблюдаемой. Так как замкнутая система должна быть асимптотически устойчивой, а преобразование Лапласа для вектора задающих воздействий в виде единичной функции 1[t ] для j - ой компоненты, может быть записано в форме: L{g}  e j 1 , где e j  (0,0,.,0,1,0,..0)T  R p , то тогда для соответствующего j - го выхода имеет s 51 место соотношение: y j ()  lim sY j ( s )  lim sWc ( s )e j s 0 s 0 1  Wc (0)e j , j  1,2,..., p , а матрица s Wc (0)  C ( A  BK ) 1 BH будет определять коэффициенты усиления замкнутой системы в установившемся режиме. Выберем матрицу H таким образом, чтобы выполнялось соотношение Wc (0)  E , то есть позиционная ошибка в установившемся режиме для каждой компоненты выхода равнялась нулю (или замкнутая система имеет астатизм первого порядка по каждой соответствующей компоненте выхода при воздействии соответствующей компоненты задающего входного воздействия, когда другие компоненты входа равны нулю). Предположим, что матрица коэффициентов усиления Wc (0) имеет размерность ( p p) . Тогда найдем, что H  [C ( A  BK ) 1 B]1 . Если rank ( H )  p , то это означает, что коэффициент усиления по какому – то выходу от соответствующего входного воздействия будет равен нулю. Если же m  p , то фактор  [C ( A  BK ) 1 B]1 будет иметь размерность ( p  m) , то есть число столбцов будет больше числа строк. В том случае, при выполнении условия rank ( H )  p можно определить матрицу H с помощью процедуры псевдо - обращения Мора – Пенроуза, которая задается соотношением: H  [C ( A  BK ) 1 B]T (C ( A  BK ) 1 B  [C ( A  BK ) 1 B]T ) 1 . В этом случае получим, что: Wc (0)  E . Восстановление фазовых координат объекта по измеряемому выходу. Построение обратной связи по состоянию осуществлялось из условия предположения, что все координаты модели состояния x(t )  R n доступны для измерения. К сожалению, для большинства реальных технических систем, данное ограничение не является выполнимым, и доступным для измерения является вектор выхода y (t )  R p . Поэтому возникает задача получения оценки x (t )  R n вектора состояния по измерениям выхода y(t ) и входа u (t )  R m . Динамическая подсистема, позволяющая, на основе измерения входа u (t ) и выхода y(t ) системы, получать оценку x (t )  R n вектора состояния, такую, что || x (t )  x(t ) || 0 при t   , называется асимптотическим наблюдателем. Если наблюдатель имеет такую же размерность, что и исходная модель состояния системы, то он называется наблюдателем полной размерности. Асимптотический наблюдатель полной размерности. Рассмотрим SISO LTI модель исходной системы: x  Ax  Bu ; y  Cx; x  R n ; u  R; y  R , у которой невозможно измерение вектора состояния x(t )  R n . Если пара матриц ( A, B) известна, то можно дублировать исходную модель, то есть построить систему вида x̂  Axˆ  Bu . Если дублирующая модель будет иметь такие же начальные условия, что и исходная система, то при любом входе u (t ) будет выполняться соотношение xˆ (t )  x(t ) для любого t  0 . Такое динамическое устройство для оценки вектора состояния исходной системы называется наблюдателем с разомкнутым контуром / 3 /. 52 Однако, такая схема построения наблюдателя имеет существенные недостатки. Например, если матрица A имеет собственные значения с положительной вещественной частью, то малейшее различие в начальных условиях x (t 0 ) и xˆ (t 0 ) будет приводить к тому, что ошибка || x (t )  x(t ) || будет расти экспоненциально во времени. Поэтому рассмотрим задачу построения асимптотического наблюдателя полного порядка с обратной связью по выходу, то основанного на ошибке между сигналом выхода исходной системы (объекта) Cx(t ) и сигналом выхода наблюдателя: Cxˆ (t ) . Обозначим вектор коэффициентов такой обратной связи, как L  R n1 , тогда схема асимптотического наблюдателя с обратной связью будет описываться уравнением: x̂  Axˆ  Bu  L( y  Cxˆ ); xˆ  R n ; u  R; y  R или x̂  [ A  LC]xˆ  Ly  Bu . Обозначим вектор ошибки оценки состояния, как e  x  xˆ . Тогда уравнения, описывающие исходную систему (объект), и наблюдатель, можно записать в виде: x  Ax  Bu; e  [ A  LC]e . Таким образом, если матрица [ A  LC ] имеет собственные значения с отрицательной вещественной частью, то все компоненты вектора ошибки оценки состояния будут стремиться к нулю, то есть: || x (t )  x(t ) || 0 при t   . Теорема. Если SISO LTI модель системы : x  Ax  Bu ; y  Cx; x  R n ; u, y  R является полностью наблюдаемой, то тогда существует подсистема асимптотического наблюдателя полной размерности с обратной связью, характеристический полином которой равен заданному n n 1  ...   0 желаемому полиному, то есть det( sE  A  LC)  o (s) , где  o ( s )  s   n 1s желаемый полином с вещественными коэффициентами. 53 Рассмотри теперь MIMO LTI модель: x  Ax  Bu ; y  Cx; x  R n ; u  R m , y  R p ; m  1 . Рассмотрим схему асимптотического наблюдателя полной размерности с обратной связью по выходу. Уравнение наблюдателя имеет вид: x̂  Axˆ  Bu  L(Cxˆ  Du  y)  ( A  LC) xˆ  ( B  LD)u  Ly , где xˆ (t ) - выход наблюдателя, L - матрица обратной связи наблюдателя, L  R n p . Обозначим ошибку оценки вектора состояния, как: e  xˆ  x . Тогда уравнение наблюдателя можно записать в форме: e  ( A  LC)e . Отсюда найдем решение: e  e ( A  LC )(t  t 0 ) e(t 0 ) . Если матрица L выбрана таким образом, что матрица ( A  LC) является гурвицевой, то тогда: lim [ xˆ (t )  x(t )]  0 . Кроме того, t  можно выбрать матрицу L таким образом, что характеристический полином: det( sE  A  LC) будет n n 1  ...   0 / 1 /. совпадать с заданным желаемым полиномом:  o ( s )  s   n 1s Теорема. Если LTI модель системы: x  Ax  Bu ; y  Cx  Du; x  R n ; u  R m , y  R p является полностью n n 1  ...   0 с вещественными наблюдаемой и, задан желаемый полином  o ( s )  s   n 1s коэффициентами, то существует асимптотический наблюдатель полной размерности: x̂  ( A  LC) xˆ  ( B  LD)u  Ly такой, что матрица [ A  LC] имеет характеристический полином: det[ sE  ( A  LC)] , который совпадает с желаемым полиномом:  o ( s )  s n   n 1s n 1  ...   0 Асимптотический наблюдатель неполной размерности. Запишем полностью наблюдаемую SISO LTI модель в виде канонической формы наблюдаемости:  b1   0 0 ...  a0      1 ...  a b    1 x  A x  B u; y  C x ; x  R n ; u, y  R , где A   , C  (0, 0, ..., 1) , B   2  .  . . ... ...  ...     0 0 ...  a  b  n 1    n Отсюда видно, что выход y(t ) совпадает с компонентой x n (t ) вектора состояния x (t ) . Следовательно, при измерении выхода, необходимо получать оценки только компонент x1 (t ),... xn 1 (t ) . Рассмотрим уравнение асимптотического наблюдателя полной размерности: x̂  A xˆ  B u  L ( y  C xˆ ); xˆ  R n ; u  R; y  R . Здесь выход наблюдателя равен: C xˆ  xˆn , а выход исходной системы (объекта) y  C x  xn  xˆ n . Выделим составляющую y  xˆn и запишем уравнение 54 ~ ~ ~ x  ( xˆ1 ,..., xˆ n 1 )T ; A x: ~ x  A~ x  Bu  a~n 1 y , где ~ наблюдателя относительно (n  1) - го вектора ~ ~ подматрица размера (n  1)  (n  1) , которая расположена в левом верхнем углу матрицы A ; B вектор (b1 , b2 ,..., bn 1 )T ; a~n 1 - вектор размера (n  1) , который совпадает с последним столбцом ~ матрицы Â . Определим следующий вектор ошибки: e~  x  ~ x , e~  R n 1 . Отсюда найдем: e~  A e~ . Рассмотрим невырожденное линейное преобразование x  Pxˆ такое, чтобы матрица A  PA P 1 имела заданный характеристический полином:  o ( s )  s n 1   n  2 s n  2  ...   0 . Для этого построим  1 0 ... 0   0     0 1 ... 0  1  матрицу P в следующем виде: P   . . ... . .  . Тогда обратная матрица равна:    0 0 ... 1   n  2   0 0 ... 0 1    1 0 ... 0  0     0 1 ... 0 1  P 1   . . ... . .  . Так как пара матриц ( A , C ) имеет каноническую форму представления,    0 0 ... 1  n  2   0 0 ... 0 1   (an 1   n  2 )  0  a0  0 0 ... 0   0    (an 1   n  2 ) 1  a1   0  1 0 ... 0  1   . . ... .  . ...  . то отсюда получим: PAP1    0 0 ... 0   n  3 (an 1   n  2 )  n  3  an  3   n  4   0 0 ... 1    n  2 ( an 1   n  2 )  n  2  an  2   n  3    0 0 ... 0  1 an 1   n  2    bn   0b1     bn 1  1b1  Соответственно, первые (n  1) компонент вектора PB имеют вид:   . При этом, ...   b   b  2 n  2 1   характеристический полином матрицы PA P 1 совпадает с заданным полиномом:  o ( s )  s n 1   n  2 s n  2  ...   0 . Теорема. Пусть SISO LTI модель системы: x  Ax  Bu ; y  Cx; x  R n ; u, y  R является полностью наблюдаемой, тогда существует асимптотический (n  1) - мерный наблюдатель, который имеет заданный характеристический полином  o ( s )  s n 1   n  2 s n  2  ...   0 . Если пара матриц ( A, C ) имеет представление в виде канонической формы наблюдаемости, то тогда выход системы представляет собой n - ый компонент вектора состояния. Такой асимптотический (n  1) - мерный наблюдатель называется наблюдателем Люенбергера. Наблюдатель Люенбергера может быть также построен и для MIMO LTI модели. 55 Теорема. Для полностью наблюдаемой MIMO LTI модели: x  Ax  Bu ; y  Cx; x  R n ; u  R m , y  R p существует асимптотический (n  p) - мерный наблюдатель, который имеет заданный характеристический полином (n  p) - го порядка, где p  rank (C ) Рассмотрим алгоритм построения наблюдателя Люенбергера (n  p) - мерной размерности. Пусть матрица C имеет ранг равный p . Выделим p линейных комбинаций y  Cx , которые определяются выходом системы. Оставшиеся (n  p) независимых компонент вектора z определим с помощью соотношения z  Lx , где матрица L  R ( n  p ) n . То есть, определим линейное  z L     x  Px . Очевидно, что в силу построения, матрица P  y C  преобразование соотношением:  z   y является обратимой (несингулярной). Тогда можно записать следующие выражения: x  P 1  z  z   PAP1   PBu . Последнее уравнение можно записать, с помощью блочной структуры  y   y и   z   Azz     y    Ayz Azy  z   Bz      u . Первое матричное уравнение имеет вид: Ayy  y   B y  z  Azz z  Azy y  Bzu , где слагаемое ( Bzyu ) y  Azy y  Bz u является измеримым, так как выход y  Cx матриц, в виде:  и вход u являются измеримыми. Тогда уравнение наблюдателя для получения оценки вектора z можно записать в виде zˆ  Azz zˆ  Azy y  Bz u; zˆ (0)  0 . В этом случае, вектор ошибки по выходу e z (t )  z (t )  zˆ (t ) будет удовлетворять уравнению: e z  Azz ez ; ez (0)  z (0)  zˆ(0) . Оценка  zˆ   . Если матрица L выбрана таким  y вектора состояния исходной системы определяется, как xˆ  P 1 образом, чтобы матрица: Azz  R ( n  p ) ( n  p ) была гурвицевой, то тогда найдем:  zˆ  e  x  xˆ  P 1    0 при t   . Структура (n  p) - мерного наблюдателя Люенбергера 0 приведена ниже на рисунке. Обычно, при построении наблюдателя не полной размерности в качестве координат выхода y используют соответствующие компоненты вектора состояния. Это приводит к упрощению расчетов. 56 Матричное уравнение TA  DT  S . Рассмотрим LTI однородную модель системы: x  Ax; y  Cx; x  R n ; y  R p . Пусть выход системы является входом наблюдателя, заданного в виде следующей динамической модели: z  Vz  LCx; z  R q , где L  R q  p , V  R qq . Будем также считать, что пара матриц (V , L) является управляемой, а матрица V является гурвицевой. Обозначим LC  S и запишем уравнения наблюдателя в виде: z  Vz  Sx . Предположим, что существует матрица линейного преобразования T  R q  n такая, что имеет место равенство: TA  VT  S , то есть z  TAx  Tx . Тогда, если матрица V не имеет общих собственных значений с матрицей A , то решение уравнения z  Vz  Sx может быть записано в форме: z (t )  Tx(t )  eVt [ z (0)  Tx(0)] . Так как матрица V является гурвицевой, то тогда: z(t )  Tx(t ) при t   для любых z (0), x(0) . Если матрица C имеет ранг равный p , то тогда решение T уравнения TA  VT  S имеет ранг равный (n  p) . Рассмотрим случай, когда матрица B  0 , то есть система имеет вид: x  Ax  Bu ; y  Cx; x  R n ; u  R m , y  R p и, пусть существует решение T уравнения TA  VT  S , где z  Tx , тогда уравнение асимптотического наблюдателя может быть записано в форме: z  Vz  Sx  TBu; z  R q . При этом, решение z имеет тот же самый вид: z (t )  Tx(t )  eVt [ z (0)  Tx(0)] . Отсюда, можно сформулировать следующее утверждение / 1 /. Теорема. Пусть задана LTI модель системы x  Ax  Bu ; y  Cx; x  R n ; u  R m и асимптотический наблюдатель, заданный уравнением: z  Vz  Sx  TBu ; z  R q , u  R m , x  R n ,где V является гурвицевой матрицей, не имеющей общих собственных значений с матрицей A ; матрица равна S  L  C и пара матриц (V , L) является полностью управляемой. Тогда выход z наблюдателя будет стремиться к вектору Tx при t   для любых начальных условий z (0), x(0) , тогда и только тогда, когда матрица T является решением матричного уравнения TA  VT  S . Пусть пара матриц pair ( A, C ); p  rank (C ) является полностью наблюдаемой. Тогда можно детализировать алгоритм построения асимптотического наблюдателя редуцированной размерности. 1. Выберем, некоторым заданием, устойчивую матрицу V  R(n  p)(n  p) , которая не имеет общих собственных значений с матрицей A . 2. Построим такую матрицу L  R ( n  p ) p , чтобы пара матриц (V , L) была управляемой. 3. Найдем решение T  R ( n  p ) n уравнения TA  VT  LC . 4. Если квадратная матрица P    порядка n является вырожденной (сингулярной), то C  T  вернуться к шагу 2. Если матрица P является несингулярной, то тогда (n  p) - мерный 1 C   y  асимптотический наблюдатель z  Vz  TBu  Ly , xˆ      позволяет получать оценку вектора T   z  состояния системы x . Матричное уравнение VX  XA  Q  0 , где матрица X является решением, часто называют расширенным матричным уравнением Ляпунова. 57 Функции среды Matlab для синтеза наблюдателей. Функция ac ker(.) позволяет рассчитывать вектор L асимптотического наблюдателя полной размерности для SISO LTI полностью наблюдаемой модели системы. Для этого используется следующий синтаксис команды: L  accer( AT , C T , p)T . Функция place(.) позволяет рассчитывать матрицу L асимптотического наблюдателя полной размерности для полностью наблюдаемой MIMO LTI модели системы. Синтаксис команды имеет вид: L  place( AT , C T , p)T , где p - вектор – строка, определяющий желаемые нули характеристического полинома матрицы ( A  LC) наблюдателя. Функция est  estim(sys, L, sensors, known) позволяет построить модель наблюдателя на основе заданной матрицы L , измеряемых компонент входа u , измеряемых компонент выхода y , не измеряемых компонент выхода z и не измеряемых возмущений w . То есть модель системы x  Ax  B1w  B2u   z D   D  . Вектора sensors и known задают индексы задается в виде:    C1       x   11 w   12 u  y   C2   D21   D22   измеряемых входов и выходов. Модель наблюдателя формируется в форме:  xˆ  Axˆ  B1w  B2u  L( y  C2 xˆ  D22u )  .  yˆ   C2   D22        xˆ   u   xˆ   E   0   Функция X  lyap(V , A, Q) позволяет решать расширенное уравнение Ляпунова. Матрицы V , A, Q должны иметь согласованные размеры, но не обязательно должны быть квадратными. Обратная связь на основе оценок состояния. Пусть задана полностью управляемая и полностью наблюдаемая LTI система: x  Ax  Bu ; y  Cx  Du; x  R n ; u  R m , y  R p и , построена матрица K обратной связи по состоянию такая, что характеристический полином матрицы замкнутой системы ( A  BK ) совпадает с заданным полиномом  d ( )  s n   n 1s n 1  ...   0 . Кроме того, пусть построен асимптотический наблюдатель полной размерности x̂  ( A  LC) xˆ  ( B  LD)u  Ly; u  g  Kx такой, что характеристический полином матрицы ( A  LC) совпадает с полиномом:  o ( s )  s n   n 1s n 1  ...   0 . Подставим в обратную связь  Kx вместо вектора x его оценку x̂ на выходе наблюдателя. Тогда уравнения замкнутой системы с наблюдателем можно записать в форме: x  Ax  BKxˆ  Bg , x̂  ( A  LC) xˆ  B( g  Kxˆ )  LCx  LDg . 58 Теорема. Пусть задана полностью управляемая и полностью наблюдаемая LTI модель системы ( A, B, C ) . Предположим, что построена матрица K обратной связи такая, что матрица ( A  BK ) имеет характеристический полином равный:  d ( )  s n   n 1s n 1  ...   0 . Предположим также, что построена матрица L асимптотического наблюдателя такая, что матрица ( A  LC) имеет n n 1  ...   0 . Обозначим через  2 n ( s ) характеристический полином равный:  o ( s )  s   n 1s характеристический полином замкнутой системы с наблюдателем, заданной уравнениями: x  Ax  BKxˆ  Bg , x̂  ( A  LC) xˆ  B( g  Kxˆ)  LCx  LDg . Тогда имеет место равенство: 2n (s)  d (s)o (s) . Доказательство. Запишем уравнения замкнутой системы с наблюдателем в следующем виде:  x   A  ̂     x   LC  BK  x   B       g . Обозначим: ~ x  x  xˆ и запишем линейное A  LC  BK  xˆ   B  LD   x E  x преобразование:  ~   M   , где M   E x  xˆ  0  , E  R n n . В новых координатах уравнения  E  x   A  BK    0 замкнутой системы примут вид:  ~    x BK  x   B       g . Отсюда видно, что A  LC  ~ x    LD  характеристический полином такой системы имеет вид: 2n (s)  d (s)o (s) . Так как построение обратной связи по состоянию и асимптотического наблюдателя осуществлялось независимо друг от друга, то не имеет значения и порядок выполнения расчетов матриц K и L . Такой подход к проектированию систем управления называется принципом разделения. 59 Рассчитаем передаточную функцию Wc (s ) замкнутой системы с наблюдателем. Имеют место следующие соотношения: Wc ( s )  C  sE  A  BK 0   BK   sE  A  LC  1  B    0  ( sE  A  BK ) 1 ( sE  A  BK ) 1 BK ( sE  A  LC ) 1  B     C ( sE  A  BK ) 1 B  P( s ) o ( s ) ,  C 0  0   d ( s) o ( s) ( sE  A  LC ) 1    P(s)  C ( sE  A  BK ) 1 B передаточная функция замкнутой системы без наблюдателя. Таким где  d (s) образом, матрица L оказывает влияние на переходной процесс в системе только при ненулевых начальных условиях x(0), ~ x (0)  0 . Так как, в передаточной функции Wc (s ) имеются сокращаемые полюса и нули, то замкнутая система с наблюдателем будет неуправляемой и/или ненаблюдаемой. Синтез линейных оптимальных систем с интегральным квадратичным функционалом. Рассмотрим следующую линейную систему: x  A(t ) x  B(t )u; y  C T (t ) x; x(t0 )  x0 . Квадратичный критерий качества определим с помощью следующего соотношения: tf J   {u T ( )u ( )  xT ( ) L( ) x( )}d  xT (t f )Qx(t f ) . Задача синтеза состоит в поиске t0 управляющей функции u (t ) , определенной на интервале t 0  t  t f , которая определена на траекториях движениях системы и минимизирует функционал J . Матрицы Q, L , обычно, полагают симметричными. Данная задача является вариационной и, часто, рассматривается в разделе теории оптимального управления. Однако, если в качестве управляющих функций рассмотреть m мерное пространство {C[t 0 , t f ]} m непрерывных функций, то можно решить указанную задачу геометрическими способами, используемыми при исследовании линейных систем /Андреев с359/. Можно выделить два вида решений приведенной задачи. В первом случае оптимальное управление u (t ) определяется как функция от матриц A, B, Q, L , начального состояния x0 ,t0 , временного параметра t и множества значений, которому должно принадлежать состояние x1 , t f . То есть, управление рассматривается для разомкнутого контура системы, полностью определяется априори заданными параметрами задачи и не зависит от текущего состояния системы. Во втором случае, управление ищется в виде линейной функции от текущего состояния ( x (t ), t ) при движении системы по оптимальной траектории. То есть ищется оптимальное управление в виде u (t )   K (t ) x(t ) . В этом случае управление представлено в виде сигнала обратной связи или, говорят, управление осуществляется по замкнутому контуру. При решении поставленных задач будем предполагать, что состояние системы, либо ее выход доступны для измерения. Рассмотрим предварительно задачу без ограничений на конечное состояние системы. Определим матричное дифференциальное уравнение, называемое уравнением Риккати, с помощью соотношения: K (t )   AT (t ) K (t )  K (t ) A(t )  K (t ) B(t ) BT (t ) K (t )  L(t ) Теорема (о существовании оптимального управления для линейной системы без ограничений на конечное состояние). Пусть заданы матрицы A(t ), B(t ), L(t ), Q(t ) , где L(t ), Q(t ) симметричные матрицы. Предположим, что на интервале t 0  t  t f существует решение K (t ) 60 дифференциального уравнения Риккати K (t )   AT (t ) K (t )  K (t ) A(t )  K (t ) B(t ) BT (t ) K (t )  L(t ) с начальным условием K (t1 )  Q . Тогда существует управление u (t ) , которое дает минимум критерию качества: tf J   {u T ( )u ( )  x T ( ) L( ) x( )}d  x T (t f )Qx (t f ) на траекториях системы t0 x  A(t ) x  B(t )u; x(t0 )  x0 . Минимальное значение J равно x T (t 0 ) K (t ) x(t 0 ) . Оптимальное управление в виде обратной связи определяется следующим соотношением: u (t )   BT (t ) K (t , Q, t f ) x(t ) . Оптимальное управление по разомкнутому контуру имеет вид: u (t )   B T (t ) K (t , Q, t f )(t , t 0 ) x0 , где (t , t 0 ) - переходная матрица для уравнения: x  [ A(t )  B(t ) BT (t ) K (t , Q, t f )]x . Решение задачи в виде обратной связи называется построением регулятора нулевого состояния. Структура замкнутой системы, при этом, будет иметь следующий вид: Пример. Пусть x (0)  1 . Необходимо найти x (t ) на интервале 0  t  T так, чтобы критерий, задаваемый T соотношением: J   {[ x ( )]2  x 2 ( )}d достигал минимума. В этом случае, уравнение движения можно t0 записать в следующем виде: x (t )  u (t ) . Тогда уравнение Риккати для этой задачи имеет вид: k(t )  k 2 (t )  1 . Решение, которое проходит через нуль при t  T , описывается функцией k (t )  th(T  t ) . Поэтому оптимальная траектория удовлетворяет уравнению: x(t )  th(T  t ) x(t ) . sh(T ) Интегрируя это уравнение, получим следующие соотношения: x(t )  ch(t )  ( ) sh(t ) , ch(T ) u(t )  th(T  t ) x(t ) . Минимальное значение квадратичной оценки качества равно: xT (0) K (0,0, T ) x(0)  x 2 (0)th(T ) . Достаточно часто, требуемое ограничение на управление, которое в явной форме не учитывается в постановке задачи оптимизации, может быть выбором соответствующей матричной функции G (t ) . Теорема. Пусть задана система x  A(t ) x  B(t )u  h(t ); y  C T (t ) x; x(t0 )  x0 ; x  R n ; u  R m ; y  R P и критерий оптимизации 61 tf вида J   {u T ( )G ( )u ( )  x T ( ) L( ) x( )}d  x T (t f )Qx (t f ) . Тогда решение задачи синтеза t0 оптимального нестационарного линейного регулятора существует и единственно, а оптимальное управление имеет вид u *  (G 1 B T Kx  матрица 1 1 T G B p) , где симметричная ( n  m ) 2 K и n - мерный вектор p определяются из уравнений K  KA  AT K  KBG 1BT K  L ; p  KBG 1BT p  AT p  2Kh при граничных условиях K (t f )  Q ; p (t f )  0 . При этом для любого t  [t 0 , t f ] справедливо соотношение tf x (t ) K (t ) x(t )  p (t ) x(t )  q (t )  x (t f )Qx (t f )   [( x * ) T L( ) x * ( )  (u * ) T G ( )u * ( )]d , где T T T t q (t ) - скалярная функция, которая определяется из уравнений 1 q  pT BG 1BT p  pT h; q(t f )  0 . 4 h(t )  0 оптимальное управление имеет вид u*  G 1BT Kx , где симметричная ( n  n ) матрица K определяется из уравнения K  KA  AT K  KBG 1BT K  L ; K (t f )  Q . Следствие. При Рассмотрим теперь случай линейной стационарной системы вида: x  Ax  Bu; y  C T x; x(t0 )  x0 ; x  R n ; u  R m ; y  R p . Теорема (о связи управляемости и идентифицируемости системы с решением алгебраического уравнения Риккати ). Если система {A, B, C} управляема и идентифицируема, то существует симметрическая положительно определенная матрица V , которая является решением алгебраического уравнения Риккати: (без доказательства). AT V  VA  VBBT V  CC T  0 Таким образом, управляемость и идентифицируемость тройки матриц . {A, B, C} является необходимым условием существования соответствующего алгебраического уравнения Риккати. Теорема (об оптимальном стационарном регуляторе). Пусть стационарная система управляема и идентифицируема. Пусть также {A, B, C} Vr - положительно определенное решение алгебраического матричного уравнения: A V  VA  VBB V  CC  0 . Тогда существует управление, которое минимизирует квадратичный критерий качества следующего вида: T  J   {u T ( )u ( )  x T ( )CC T x( )}d для системы T T x  Ax  Bu; y  C T x; x(t0 )  x0 ; . Причем t0 минимальное значение имеет вид: J равно x0T Vr x0 , а оптимальное u(t )   BT Vr x(t ) , и в разомкнутой форме: управление в виде обратной связи T V ]t r u (t )   B T Vr e[ A BB . 62 Доказательство. Сделаем замену переменной: v(t )  u (t )  B Vr x(t ) , называемую T преобразованием Риккати. Уравнения движения примут вид: x (t )  [ A  BB Vr ] x(t )  Bv (t ) , и T функционал J будет определяться соотношением:  J   [v( )  B T Vr x( )]T {[v( )  B T Vr , x( )]  x T ( )CC T x( )}d . Преобразуем подынтегральное t0 выражение, применив равенство: AT Vr  Vr A  Vr BB T Vr  CC T . Получим, после всех преобразований, следующее соотношение:  J   {[v T ( )v( )  [( A  BB T Vr ) x( )  Bv ( )]T Vr x( )  x T ( )Vr [( A  BB T Vr ) x( )  Bv ( )]}d t0    или J   v ( )v ( )d  x (t )Vr x (t ) |t0 . Очевидно, что минимум значения J достигается при T T t0 v(t )  u (t )  B T Vr x(t )  0 . Отсюда следует доказательство утверждения теоремы. Структурная схема стационарного оптимального регулятора приведена на рисунке. Имеет место также и более общий результат. Пусть линейная стационарная система описывается уравнением: x  Ax  Bu , а критерий оптимальности имеет вид:  J   {u T ( )Gu ( )  x T ( )Qx( )}d (*) Теорема (об оптимальном стационарном регуляторе). Задача (*) имеет решение тогда и только тогда, когда стационарная система {A, B} стабилизируема, а оптимальное управление имеет вид u(t )  G 1BT Vr x(t ) , где алгебраического матричного уравнения: - положительно определенное решение Vr AT V  VA  VBG 1BT V  Q . При этом для любого t t  0 справедливо равенство x (t )Vr x(t )   {u T ( )Gu ( )  xT ( )Qx( )}d . T