Математическое и компьютерное моделирование
 Скачать 3.02 Mb.
|
|
2.3. Дискретизация непрерывных моделей 2.3.1. Решение уравнений состояния Пусть исследуемая система линейна и не содержит нестационарных элементов, т.е. уравнения динамики имеют вид ), ( ) ( ) ( t Bu t Ax dt t dx + = (2.72) ), ( ) ( ) ( t Du t Cx t y + = (2.73) где x(t) ∈R n ; u(t) ∈R m ; y(t) ∈R l . Очевидно, что если функция x(t) определена, то выход y(t) вычисляется непосредственно из (2.73). Рассмотрим решение задачи Коши, т.е. определим вектор состояния x(t) при известном начальном значении х(0) и заданном входном процессе и(t). Решение однородного уравнения ) ( ) ( t Ax dt t dx = (2.74) для всех действительных t 0 , при выполнении равенства , ) , ( Ф ), , ( Ф ) , ( Ф 0 0 0 0 E t t t t A dt t t d = = (2.75) 42 где x(t) ∈R n ; Ф(t, t 0 ) = X(t)Х -1 (t 0 ) – переходная матрица состояния; X(t) – фундаментальная матрица-функция,можно представить в виде ). ( ) , ( Ф ) ( 0 0 t x t t t x = (2.76) Для того, чтобы воспользоваться выражением (2.76), необходимо уметь вычислять матрицу Ф(t, t 0 ). Известно, что поскольку уравнение (2.74) стационарное, то матрица Ф(t, t 0 ) зависит от одного аргумента θ и совпадает с матричной экспонентой Ф(t, t 0 ) = е А θ , θ = t – t 0 , которую мож- но определить в виде ряда ! ) ( ! ) ( 2 ) ( 1 2 ∑ ∞ = + = + + + + + + = k k k A k A E k A A A E e θ θ θ θ θ (2.77) Таким образом, решение уравнения (2.74) определяется формулой , ) ( 0 x e t x At = (2.78) где х 0 = х(t 0 =0). Выражение (2.78) имеет наиболее простой вид тогда, ко- гда матрица А = diag{s 1 , s 2 , …, s n }, поскольку матрица e At также будет диагональной { } t s t s t s At n e e e e diag 2 1 = Решение неоднородного уравнения (2.72) можно представить в виде ), ( ) ( ) ( t x t x t x част общ + = где x общ (t) – переходная составляющая или общее решение однородного уравнения (2.74) при заданных начальных условиях; x част (t) – вынужден- ная составляющая или частное решение уравнения (2.72) при нулевых начальных условиях. Как известно, x общ (t) имеет вид ∫ = t t общ d Bu t t x 0 , ) ( ) , ( Ф ) ( ϑ ϑ ϑ а x част (t) представляет собой выражение (2.76); следовательно, для x(t) можно записать следующую формулу Коши , ) ( ) , ( Ф ) , ( Ф ) ( 0 0 0 ∫ + = t t d Bu t x t t t x ϑ ϑ ϑ (2.79) а для стационарных систем в силу Ф(t, t 0 ) = е А θ , θ = t – t 0 выражение (2.79) представить в виде 43 ( ) ( ) , ) ( ) ( 0 0 0 ∫ − − + = t t t A t t A d Bu e x e t x ϑ ϑ ϑ (2.80) часто называемом переходным уравнением состояния. Свойства переходной матрицы. Перечислим основные свойства мат- рицы Ф(t, t 0 ). 1. ∀t 0 имеет место Ф(t 0 , t 0 ) = E. 2. ∀t 0 , t 1 , t выполнено правило композиции Ф(t, t 0 ) = Ф(t, t 1 )Ф(t 1 , t 0 ). 3. det Ф(t, t 0 ) ≠ 0 ∀t 0 , t. 4. Ф(t, t 0 ) = X(t)Х -1 (t 0 ) , где X(t) – любая фундаментальная матрица. 5. (Ф(t, t 0 )) -1 = Ф(t 0 , t) ∀t 0 , t. 6. Справедливо уравнение ) , ( ), , ( ) ( ) , ( 0 0 0 0 E t t Ф t t Ф t A dt t t dФ = = 7. Матрица Ф T (t, t 0 ) удовлетворяет сопряженному уравнению ) , ( ), , ( ) ( ) , ( 0 0 0 0 E t t Ф t t Ф t A dt t t dФ T T T = − = 8. Если det T ≠ 0, Ф T (t, t 0 ) = Т -1 Ф T (t, t 0 )Т, где Ф (t, t 0 ) удовлетворя- ет (2.75), в котором вместо матрицы А подставлена подобная ей мат- рица A(t) = ТA(t)Т -1 2.3.2. Дискретные модели непрерывных систем Постановка задачи дискретизации. Для динамической системы, не- прерывная модель которой имеет вид , ), ( ) ( ) ( ), ( ) ( ) ( R t t Du t Cx t y t Bu t Ax dt t dx ∈ + = + = (2.81) требуется получить эквивалентную систему разностных уравнений или ее дискретную модель вида ,... 1 , 0 ], [ ] [ ] [ ], [ ] [ ] 1 [ = + = + = + k k u D k x C k y k Qu k Px k x (2.82) Здесь под эквивалентностью систем понимается совпадение их ре- 44 акций на одно и то же входное воздействие при соответствующих на- чальных условиях. Более того, при u[k] = u(t k ), где t k = kT 0 , T 0 = const – шаг дискретизации (период квантования), выполнено у[k] = у(t k ) – решения уравнений (2.81) и (2.82) совпадают при t k = kT 0 Формулы перехода к разностным уравнениям. Рассмотрим задачу определения матриц D C Q P , , , в (2.82) по известным матрицам А, B, C, D в (2.81), исходя из сформулированного выше требования эквивалентно- сти систем по отношению к входному воздействию u(t). Для простоты изложения ограничимся рассмотрением кусочно-постоянного процесса вида , ,... 2 , 1 , 0 , , при ) ( ) ( 0 1 = = < ≤ = + k kT t t t t t u t u k k k k достаточно распространенного на практике, который формируется так называемым фиксаторам или экстраполяторам нулевого порядка. Из- вестно, что решение этой задачи с использованием аппарата передаточ- ных функций и z-преобразования (дискретного преобразования Лапласа) позволяет передаточную функцию дискретной модели W D (z) записать следующим образом: ( ) , ) ( 1 ) ( 1 − = − s s W Z z z W (2.83) где Z означает операцию z-преобразования переходной функции исход- ной непрерывной системы. Для решения аналогичной задачи в рамках метода пространства со- стояний воспользуемся формулой Коши (2.80) и проинтегрируем первое уравнение из (2.81) на интервале [t k , t k+1 ], полагая на нем u(t)= u(t k ). При х 0 = х(t k ) получим ( ) ( ) ( ) ). ( ) ( ) ( ) ( ) ( 1 1 0 1 1 1 1 k t t t A k AT t t t A k t t A k t Bu d e t x e d Bu e t x e t x k k k k k k k k + = = + = ∫ ∫ + + + + + − − − + υ υ υ υ υ (2.84) Если ввести новую переменную θ = t k+1 – υ , то, вычисляя интеграл в круглых скобках из соотношения (2.84), при t k+1 – t k = Т 0 получим , ( ) ∫ ∫ − = = − − + + 0 0 1 1 0 1 ) ( T AT A t t t A E e A d e d e k k k θ υ θ υ 45 предполагая, что матрица А невырожденная. Следовательно, выражение (2.84) можно переписать в виде , 0 det ), ( ) ( ) ( ) ( 0 0 1 1 ≠ − + = − + A t Bu E e A t x e t x k AT k AT k (2.85) а уравнение выхода, согласно второму уравнению из (2.81), представить в виде ). ( ) ( ) ( k k k t Du t Cx t y + = (2.86) Сопоставляя первое уравнение из (2.82) с уравнением (2.85) и вто- рое уравнение из (2.82) с уравнением (2.86), получим очевидные соотно- шения , , ) ( , 1 0 D D C C B E P A Q e P AT = = − = = − (2.87) Когда выполнен переход от (2.81) к (2.82), то с учетом соотношений (2.87) можно записать уравнения , , 1 , 0 ], [ ] [ ] [ , 0 det ], [ ) ( ] [ ] 1 [ 0 0 1 = + = ≠ − + = + − k k Du k Cx k y A k Bu E e A k x e k x AT AT (2.88) которые с помощью z-преобразования и преобразований, аналогичных соотношениям (2.67), (2.68), могут быть представлены следующим обра- зом: ( ) ( ) ( ) ( ) 0 det ), ( ) ( ) ( ) ( ) ( 1 1 1 0 0 ≠ + − = = + − − = = = − − − A z u D Q P zE C z u D B E e A e zE C z u z W z y AT AT (2.89) Заметим, что передаточные функции в выражениях (2.83) и (2.89) идентичны, но получены различными способами. 2.3.3. Вычисление матричной экспоненты Методы вычисления матричной экспоненты e At обычно подразделя- ют на точные и приближенные. Точные позволяют получить выражение для матричной экспоненты через скалярные аналитические функции, а приближенные – аппроксимировать ее с некоторой алгоритмической ошибкой, зависящей от способа аппроксимации и параметров алгоритма. 46 Точные методы. Рассмотрим аналитические выражения для матрич- ной функции e At в следующих частных случаях: 1. Матрица А диагональная с вещественными собственными значе- ниями. Пусть А = diag{s 1 , s 2 , …, s n }, Im s i = 0, i = 1, …, n. Непосредствен- ным вычислением суммы ряда (2.77) получаем { } t s t s t s At n e e e e diag 2 1 = , где e – скалярные экспоненты. t s i 2. Матрица А блочно-диагональная с мнимыми собственными зна- чениями. Пусть А имеет собственные числа: s 1,2 = α ± j β , j 2 = –1, тогда ее можно представить в виде − + = − = 0 0 β β α α β β α E A , а матричную экспоненту описать выражением 0 0 t Et At e e e − = β β α Применяя формулу (2.77), можно показать, что это соотношение имеет вид cos sin sin cos − = t t t t e e t At β β β β α 3. Матрица А имеет кратные вещественные собственные значения. Пусть s i = 0, i = 1,2,3, т.е. матрица А нильпотентная , 0 0 0 1 0 0 0 1 0 = A вычисляя степени которой получаем, что 0 , 0 0 0 0 0 0 1 0 0 4 3 2 = = = = A A A Следовательно, ряд (2.77) точно выражается первыми тремя слагае- мыми в виде 47 1 0 0 1 0 5 0 1 2 = t t t e At В случае отличных от нуля кратных вещественных собственных зна- чений s i = α , i = 1,2,3 аналогично предыдущему получаем 1 0 0 1 0 5 0 1 2 = t t t e e t At α Если матрица А имеет произвольный вид, то всегда существует не- вырожденное преобразование с матрицей Т такое, что подобная ей мат- рица A = ТАТ -1 – жорданова. Тогда по свойству 8 переходной матрицы (см. 2.3.1) получаем T e T e t A At 1 − = Аналитические формулы для матричной экспоненты могут быть по- лучены также на основе преобразования Лапласа. Этот способ основан на том, что резольвента R(s) постоянной матрицы А является изображением по Лапласу ее матричной экспоненты, т.е. ( ) ( ) , 1 − − = A sE e L At причем элементы переходной матрицы можно найти с помощью таблиц обратного преобразования Лапласа. Приближенные методы. Эти методы основаны на различных ап- проксимациях ряда (2.77) выражениями, содержащими конечное число слагаемых. Достаточно очевидной является аппроксимация Тейлора порядка k, согласно которой ряд (2.77) приближенно заменяется следующей конеч- ной суммой: ! ) ( ! ) ( 2 ) ( 1 2 ∑ = + ≡ + + + + ≈ k i i k A i A E k A A A E e θ θ θ θ θ (2.90) В частности, при k = 1 получаем линейное приближение вида , θ θ A E e A + ≈ (2.91) которое будем называть аппроксимацией Эйлера. Аппроксимация (2.90) не является лучшей и во многих отношениях 48 уступает более общей аппроксимации Паде. При такой аппроксимации экспонента е х представляется рациональной функцией , ) ( ) ( x G x F e x µν µν ≈ (2.92) ( ) ( ) ( )( ) ( ) ( )( ) ( ) , ! 1 1 1 2 1 ! 2 1 1 ! 1 1 ) ( 2 µ µν µ ν ν µ ν µ µ µ ν µ ν µ µ µ ν µ µ x x x x F + ⋅⋅ ⋅ − + + ⋅ ⋅⋅ ⋅ − + + − + + − + + + = (2.93) ( ) ( ) ( )( ) ( ) ( )( ) ( ) ! 1 1 1 2 1 ) 1 ( ! 2 1 1 ! 1 1 ) ( 2 ν ν µν ν ν ν µ ν µ ν ν ν µ ν µ ν ν ν µ ν x x x x G + ⋅⋅ ⋅ − + + ⋅ ⋅⋅ ⋅ − − + + − + + − + + − = (2.94) Соответственно, для матричного аргумента x = A θ запишем ), ( ) ( 1 θ θ µν µν θ A G A F e A − ≈ (2.95) где F µν (A θ ), G µν (A θ ) – матричные многочлены вида (2.93), (2.94). В даль- нейшем (2.95) будем называть аппроксимацией Паде ( µ , ν ). Приведем некоторые частные случаи (2.95). Во-первых, аппроксимация Тейлора (2.90) – это частный случай (2.95) при ν = 0. Следовательно, выражение (2.91) совпадает с аппрокси- мацией Паде (1, 0). Аппроксимация Паде (0, 1) имеет вид ( 1 − − ≈ θ θ A E e A ) ) (2.96) и далее будет называться неявным методом Эйлера. Аппроксимация Паде (1, 1) соответствует методу Тастина и опреде- ляется формулой ( )( 0.5 0.5 1 − − + ≈ θ θ θ A E A E e A (2.97) Формула Паде (2, 2) приводит к выражению ( )( ) ) ( 6 12 ) ( 6 12 1 2 2 − + − + + ≈ θ θ θ θ θ A A E A A E e A (2.98) Основным преимуществом аппроксимаций Паде является их более высокая точность, а недостатком – необходимость обращения матрицы 49 G µν (A θ ) и связанная с этим проблема ее вырожденности. 2.3.4. Вычисление матрицы Q Как уже отмечалось, формула (2.87) для вычисления матрицы Q применима, если det A ≠ 0. Трудностей, связанных с вычислением Q при вырожденной матрице А, можно избежать, если при формальной подстановке выражения 0 AT e P = , полученного из аппроксимаций Тейлора (2.90) или Паде (2.95), в соотношении (2.87) выполнить "сокращение" матрицы А. Тогда в выра- жение для Q матрица А -1 входить не будет. В частности, аппроксимация по методу Эйлера (2.87), (2.91) Р = Е + АТ 0 приводит к формуле Q = ВТ 0 , а аппроксимация Паде (1, 1) (2.87), (2.97) – к формуле Q = (Е – 0.5АТ 0 ) -1 ВТ 0 Иной способ заключается в расширении пространства состояний ис- ходной системы (2.81). Входной процесс и(t) при t k ≤ t < t k+1 рассматрива- ется как решение некоторого однородного дифференциального уравне- ния. Тогда расширенная система тоже однородная и в вычислении по (2.87) нет необходимости. Искомые матрицы Р и Q получаются как под- матрицы расширенной "матричной" экспоненты. Покажем это на приме- ре ступенчатого входного процесса и(t) = и(t k ) при t k ≤ t < t k+1 . Для ука- занного промежутка времени уравнение состояния (2.81) запишем в виде ]. [ ) ( , 0 ) ( , ], [ ) ( ), ( ) ( ) ( 1 k u t u dt t du t t t k x t x t Bu t Ax dt t dx k k k k = = < ≤ = + = + (2.99) Введем в рассмотрение расширенный (n + m)-мерный вектор со- стояния )} ( ,) ( { col ) ( t u t x t x = и (n + m)х(n + m) матрицу 0 0 x x = m m m n B A A Уравнение (2.99) запишем в виде ]}, [ ], [ { col ) ( ), ( ) ( 1 + < ≤ = = k k k t t t k u k x t x t x A dt t x d (2.100) Соответствующая дискретная модель, аналог (2.82), принимает вид ], [ ] 1 [ 0 T A e P k x P k x = = + (2.101) 50 Учитывая структуру матрицы A и формулу (2.77), для матрицы P непосредственно находим 0 ' ' = E Q P P С учетом этого из (2.101) получаем: ]. [ ' ] [ ' ] 1 [ k u Q k x P k x + = + (2.102) Сравнивая (2.102) с (2.82), имеем очевидный факт – матрицы Р и Q в (2.82) совпадают с матрицами Р' и Q'. Другими словами, они могут быть получены как соответствующие подматрицы матрицы P вида (2.101). 2.3.5. Вычисление передаточной функции дискретной модели Уже отмечалось, что уравнение (2.89) описывает дискретную систе- му с помощью передаточной функции, вычисленной по разностному уравнению (2.82), исходя из преобразования уравнений состояния непре- рывной системы (2.81). Однако можно получить и приближенное решение задачи, при кото- ром искомая функция W(z) определяется непосредственной заменой ар- гумента s в W(s). Эти формулы основаны на "линейных" аппроксимациях Паде ( µ , ν ), в которых значения µ и ν не превышают единицы. Для формулы аппроксимации Эйлера (2.91) в (2.89) следует подста- вить Р = Е + АТ 0 и, как показано в п. 2.3.4, учесть Q = ВТ 0 . Таким обра- зом, получаем ( ) ( ) 1 ) 1 ( ) ( 1 0 0 1 0 1 B A T z C BT AT E z C Q P zE C z W − − − − − = = − − = − = (2.103) Теперь при сравнении выражения (2.103) с формулой W(s) = C(sE – A) -1 B становится вполне очевидным, что W(z) можно получить прибли- женно из W(s) с помощью следующей замены аргумента: 0 1 ) ( ) ( T z s s W z W − = = (2.104) Для формулы неявного метода Эйлера (2.96), аналогичным образом получаем 51 0 1 ) ( 1 ) ( zT z s s W z z W − = = (2.105) Для формулы с аппроксимацией Паде (1, 1) после некоторых преоб- разований получаем подстановку метода Тастина 1 1 2 0 ) ( 1 2 ) ( + − = + = z z T s s W z z W (2.106) Точность приведенных методов зависит от соотношения между ин- тервалом Т 0 и наименьшей постоянной времени непрерывной системы W(s). 2.3.6. Устойчивость дискретных моделей Процедура построения дискретных моделей непрерывных систем всегда сопровождается требованием о сохранении свойства устойчиво- сти: устойчивая непрерывная система должна приводить к устойчивой дискретной одели, а в случае неустойчивой исходной системы ее дис- кретная модель тоже должна получиться неустойчивой. м Для точных методов перехода это требование выполнено, а для при- ближенных не всегда. Точный подход. Как известно, асимптотическая устойчивость ли- нейных непрерывных систем имеет место, если корни характеристиче- ского многочлена (собственные числа) матрицы А в (2.81) имеют отрица- тельные вещественные части, т.е. Re s i < 0 при det(s i E – A) = 0, i =1,…, n. При этом дискретная система (2.82) будет асимптотически устойчи- вой, если корни ее характеристического многочлена удовлетворяют усло- вию z i < 1 при det(z i E – P) = 0, i =1,…, n. Известно следующее утверждение: согласно точной формуле (2.87) 0 AT e P = , следовательно, во-первых, собственные числа z i матрицы Р оп- ределяются соотношением , i =1,…, n, во-вторых, при всех Т 0 T s i i e z = 0 > 0 выполнено {Re s i < 0 ⇔ z i < 1}, в-третьих, свойства устойчивости сис- тем (2.81) и (2.82) эквивалентны. Метод Эйлера. Согласно этому методу по (2.91) получаем равенство 52 Р = Е + АТ 0 . Следовательно, z i = 1+ s i Т 0 , i =1,…, n. Проверка условия устойчивости z i < 1 приводит к следующим неравенствам: ( ) ,..., 1 , Im , Re , 1 1 0 0 2 2 n i s T s T i i i i i i = = = < + + β α β α (2.107) Условие (2.107) означает, что значения корней характеристического многочлена системы, умноженные на интервал квантования, должны на- ходиться на комплексной плоскости внутри круга единичного радиуса с центром в точке (–1, j0). Это эквивалентно неравенству Re min 2 2 0 i i i s s T = (2.108) Область применения явного метода Эйлера существенно ограничи- вается малыми (относительно модулей собственных чисел системы) зна- чениями Т 0 Неявный метод Эйлера. Согласно формуле (2.96) Р = (Е – А θ ) -1 име- ем , ,..., 1 , 1 1 0 n i T s z i i = − = следовательно, условие z i < 1 приводит к неравенствам ( ) ,..., 1 , Im , Re , 1 1 0 0 2 2 n i s T s T i i i i i i = = = > + − β α β α (2.109) Условие (2.109) означает, что корни характеристического многочле- на системы, умноженные на интервал квантования, должны находиться на комплексной плоскости вне окружности единичного радиуса с цен- тром в точке (1, j0). Из этого следует, что при такой аппроксимации для любой устойчи- вой непрерывной системы будет получена устойчивая дискретная модель при всех (а не только малых) Т 0 > 0. Следует отметить, что свойства устойчивости непрерывных и дис- кретных моделей не будут эквивалентными: устойчивые дискретные мо- дели могут получиться и для неустойчивых исходных непрерывных сис- тем. Точность аппроксимации по этому методу невелика. Аппроксимация Паде при ν = µ . Можно показать, что для аппрокси- маций Паде ( ν , µ ) при ν = µ , Т 0 > 0 справедливо {Re s i < 0 ⇔ z i < 1}. Таким образом, при их использовании устойчивость системы (2.81) экви- валентна устойчивости системы (2.82). В частности, это свойство выпол- 53 нено для рассмотренных выше аппроксимаций (2.97) и (2.98). Устойчивость методов численного интегрирования. Можно устано- вить тесную связь между рассмотренными методами построения дис- кретных моделей непрерывных системи методами численного решения задач Коши. Пусть требуется проинтегрировать уравнение 0 , ) 0 ( ), , ( ) ( 0 ≥ = = t x x t x f dt t dx (2.110) Известно, что его решение может быть получено в виде некоторого рекуррентного соотношения x k+1 = ϕ (x k , k), где k = 0, 1, 2, … – номер шага (итерации), а значения x k соответствуют значениям искомой функции x(t) в дискретные моменты времени t k = kh, где h > 0 – шаг интегрирования. Вид функции ϕ (x k , k) определяется по искомой функции f(x, t) согласно выбранному методу численного интегрирования. Например, используя метод Эйлера, получаем хорошо известную формулу ,... 2 , 1 , 0 , , ) , ( 1 = ⋅ = + = + k h k t h t x f x x k k k k k (2.111) Более подробно рассмотрим случай, когда функция f(x, t) линейна по x и уравнение (2.110) может быть представлено в виде ), ( ) ( t Ax dt t dx φ + = (2.112) где А – nxn-матрица. Формула метода Эйлера (2.111) приводит к разност- ному уравнению ( ) ) ( 1 h t x Ah E x k k k φ + + = + (2.113) Обратим внимание на то, что уравнение (2.113) с точностью до обо- значений совпадает с уравнением (2.82), в котором матрицы P, Q получе- ны на основе приближенной формулы (2.91) для матричной экспоненты. Устойчивость полученной численной процедуры определяется при- веденными в п. 2.3.1 свойствами данной аппроксимации: так, для А и Т 0 = h должно выполняться неравенство (2.107). 54 Следовательно, если собственные числа системы велики по модулю (что соответствует малым по величине постоянным времени), то для по- лучения устойчивого решения необходимо выбирать шаг интегрирования h достаточно малым. Это приводит к значительным затратам машинного времени, а также надо иметь в виду, что уменьшение величины h вызыва- ет увеличение инструментальных ошибок, связанных с конечностью раз- рядной сетки ЦВМ. Подобным свойством обладают все так называемые явные методы численного интегрирования, которые для линейных стационарных систем сводятся к аппроксимации Тейлора (2.90). Отметим, что недостатки явных методов заметно проявляют себя при решении жестких систем, у которых собственные значения отлича- ются друг от друга по модулю на несколько порядков. 55 |