Моделирование разностных уравнений

Цель работы:ознакомиться с методикой моделирования разностных уравнений с помощью пакета MATLAB и SIMULINK.

 

ТЕОРЕТИЧЕСКИЕ СВЕДЕНИЯ

 

Линейное однородное разностное уравнение второго порядка с постоянными коэффициентами имеет вид:

.                           (1)

 

Чтобы подчеркнуть дискретный характер изменения времени, это уравнение часто записывают в форме

.

 

Разностные уравнения третьего и более высоких порядков записываются аналогично.

Известны два основных метода решения линейных разностных уравнений – с помощью характеристического полинома и с использованием z -преобразования, аналогичного преобразованию Лапласа.

 

1.1. Решение разностных уравнений с помощью характеристического полинома

Будем искать решение однородного разностного уравнения в виде x (t)= z^t, где z – некоторое число. Подставляя x (t) в разностное уравнение (1) при   f (t)=0  и сокращая на z^t, получаем характеристическое уравнение

 

Если его корни z ₁, z₂  вещественные и различные, то общее решение имеет вид:

 

.

 

Если   z ₁= z₂,, то в решении появляется линейный множитель

 

 

В случае пары комплексно-сопряженных корней   z _1,2= a ± i b  решение может быть записано в вещественной форме

 

, .

 

Здесь r – модуль комплексного числа   z ₁,  а   j – его аргумент.

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

Пример 1.    Решим разностное уравнение x (t +2) –5 x (t +1) + 6 x (t) = 0.

Характеристическое уравнение имеет вид z ² –5 z + 6 = 0.

Его корни вещественные и различные: z ₁ = 3, z ₂ = 2.

Общее решение: x (t) = c ₁ 3 ^t + c ₂ 2 ^t.

Пример 2.    Решим разностное уравнение x (t +2) + 2 x (t +1) + 4 x (t) = 0.

Характеристическое уравнение: z ² + 2 z + 4 = 0. У него комплексные корни: . Их положение на комплексной плоскости   z = α + i b  показано на рис. 1.

Модуль и аргумент корней можно найти непосредственно на рис.1: .

Общее решение:

Произвольные постоянные c_i находят, задавая начальные условия. Пусть, например, в примере 2 заданы начальные условия   x (0)=2; x (1)= –4.  Записываем общее решение для   t =0  и   t =1:

Отсюда находим   с₁ =2,   с ₂=0. Следовательно, решение имеет вид

Общее решение неоднородного разностного уравнения представляет со­бой сумму общего решения соответствующего однородного уравнения и любого частного решения неоднородного уравнения. Частное решение ищут в том же виде, что и правая часть, т.е. функция   f (t)  в уравнении (1):

– если f (t) – постоянная, то в виде константы;

– если f (t) – экспонента, то в виде экспоненты с тем же показателем;

– если f (t) =sin kt или cos kt, то в виде   c ₁sin kt + c ₂cos kt.

Коэффициенты   с ₁  и   с ₂  находят, подставляя частные решения в  разностное уравнение и приравнивая одноименные функции справа и слева.

Пример 3. Дано неоднородное разностное уравнение второго порядка

 

 

Находим корни характеристического полинома

 

Частное решение ищем в виде x _част= c. Подставляя его в исходное уравнение, находим, что x _част=2.

Общее решение неоднородного уравнения получаем как сумму частного решения и общего решения соответствующего однородного уравнения

 

 

Коэффициенты с ₁, с ₂ находим из уравнений , откуда .

 

1.2. Решение разностных уравнений с помощью z-преобразования

При описании дискретных систем и решении разностных уравнений широко применяется аппарат z -преобразования – это дискретный аналог преобразования Лапласа. Например, умножение изображения F (p) на оператор Лапласа p соответствует дифференцированию непрерывной функции f (t). Умножение изображения F (z) на оператор    z соответствует сдвигу функции f (t) (которая может быть непрерывной, дискретной или решетчатой) на один такт.

Таким образом, если операторы р и 1/ р – это операторы дифференцирования и интегрирования, то операторы z и z ^-1 – это операторы сдвига влево и вправо. С инженерной точки зрения оператор z ^-1 представляет собой элемент задержки.

Существуют также определенные параллели между изображениями функций F (p) и F (z). Например, изображению по Лапласу   F (p) = 1 соответствует дельта-функция   f (t) = d (t), а изображению F (z) = 1, соответствует единичный импульс. В том и другом случае оригиналом является  элементарное импульсное воздействие.

Краткая таблица z -преобразований других функций приведена в Приложении.

Пусть дано разностное уравнение n -го порядка

 

y (t+n)+ a_{n -1} y (t+n –1)+... + a ₀ y (t)= f (t)                                               (2)

 

с начальными условиями y (0)= y ₀; y (1)= y ₁;...; y (n –1)= y_{n - 1}. 

Алгоритм его решения с помощью z -преобразования содержит 4 шага.

Шаг 1. Применить z -преобразование к уравнению (2), заменяя f (t) на F (z), y (t) на Y (z); y (t +1) на z (Y (z)– y ₀) и т.д.

Шаг 2. Из полученного алгебраического уравнения выразить Y (z).

Шаг 3. Выполнить разложение Y (z) на простые дроби.

Шаг 4. Пользуясь таблицей, выполнить обратное z -преобразование.

Результатом будет искомое решение разностного уравнения.

Пример 4. Требуется решить разностное уравнение второго порядка

 

y_{n +2} –5 y_{n +1}+6 y_n = u_n

 

с нулевыми начальными значениями y ₀, y ₁и входным сигналом u_n =1.

Шаг 1. Применяем к нему z -преобразование

 

Шаг 2. Выражаем Y (z) и подставляем :

 

Шаг 3. Представляем правую часть в виде суммы простых дробей с переменной   z в числителе:

.

 

В отчете разложение делается методом неопределенных коэффициентов. В пакете MATLAB это можно сделать с помощью команды residue.

Шаг 4. С помощью таблицы z -преобразований или  команды i ztrans тулбокса SYMBOLIC пакета MATLAB находим оригиналы каждого из слагаемых и складываем их:   

y_n =0,5–2 ⁿ +0,5×3 ⁿ.

 

1.3. Система линейных разностных уравнений

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

 

X(t +1)=AX(t),      X(0)=X₀,                                            (3)

 

где  – вектор переменных, А – квадратная матрица постоянных коэффициентов,
X₀ – вектор начальных условий.

Решение этой системы может быть записано в компактной степенной форме

X(t)=A ^t X₀.                                                                      (4)

При работе в пакете MATLAB этот способ удобнее всего.

Пример 5. Найдем решение системы разностных уравнений второго порядка

 

x (t +1) =  2 x (t),  

         y (t +1) = –2 x (t) +2 y (t),

если .

Матричное описание задачи имеет вид

 X(t +1)=AX(t),      ,       .

 

Для получения решения воспользуемся степенной формулой. Последовательно возводя матрицу А в степени 2, 3,..., n (проделайте это!), находим

 

A ⁿ = 2 ⁿ , X(t)=A ^t X₀=

 

Переходя к скалярной форме записи, получаем окончательный ответ   x = 2 ^t,   y = (20– t) 2 ^t .

Другой путь решения этой задачи связан с переходом от системы двух разностных уравнений первого порядка к одному уравнению второго порядка относительно y:

y (t +2) – 4 y (t +1) + 4 y (t) = 0.

 

Оно получается, если выразить из второго уравнения   x (t),   x (t +1) и подставить в первое.

 

1.4. Моделирование дискретных систем в MATLAB и SIMULINK

В MATLAB при моделировании линейных системы с дискретным временем используются разностные уравнения, модели в пространстве состояний и дискретные передаточные функции. В ядре MATLAB имеется команда filter, позволяющая рассчитать выходной сигнал фильтра по известному входному сигналу, заданному массивом своих значений. Ее входными аргументами являются векторы а и b коэффициентов  разностного уравнения, а также массив значений сигнала и: y=filter(b,a,u). С помощью четвертого входного аргумента можно задавать начальные условия фильтра.

Дискретная система в пространстве состояний, как и непрерывная система, задается четверкой матриц:

 

                                 (5)

 

Здесь – вектор состояния; u_k, y_k – входной и выходной сигналы, A, B, C, D – постоянные матрицы.

Дискретной передаточной функцией линейной системы (5) называется отношение z -преобразования выходного сигнала к z -преобразованию входного сигнала при нулевых начальных условиях Q(z) = Y(z)/U(z). Для систем с одним входом и одним выходом она представляет собой отношение двух полиномов

 

 

Дискретная передаточная функция системы, заданной  описанием в пространстве состояний (5), может быть найдена по формуле Q(z)=C(z E-A)^–1B+D.

Дискретные модели в MATLAB задаются при помощи тех же конструкторов ss, tf и zpk, что и непрерывные с той разницей, что последним параметром задается частота дискретизации ts(sampling time).

Пример 6. Рассмотрим дискретную систему, описываемую уравнениями

 

при шаге дискретизации, равном 1.

Создадим соответствующую ss-модель с дискретным временем:

 

>> a=[-1 -1/2;1/2 -1]; b=[1;0]; c=[0 1]; d=0; ts=1.0; s1=ss(a,b,c,d, ts).

 

На дисплей будут выведены матрицы:

 

a = -1     -0.5   0.5 -1

b = 1       0

c = 0 1

d =  0

 

и сообщение: Sampling time: 1.0 Discrete-time model.

 

Аналогично создаются tf и zpk-модели с дискретным временем. Обращение к полям дискретной модели осуществляется так же, как и у непрерывной, например для вывода матрицы c надо набрать код s1.c. Время дискретизации хранится в поле ts.

Для моделирования дискретных систем используются те же команды (step, impulse, initial, lsim), что и в непрерывном случае. При получении весовой функции дискретных систем вместо реакции на дельта-функцию рассматривают реакцию на единичный импульс вида . Для этого используется та же команда impulse. Реакцию на единичный скачок (переходную функцию) получают при помощи команды step.

В тулбоксе SYMBOLIC имеются команды ztrans и i ztrans для выполнения прямого и обратного z –преобразований.

Пример. Дискретная система задана уравнениями

 

 

Найдем двумя способами (символьным и численным) передаточную функцию от входа u до выхода у.

В данном случае имеем следующие матрицы описания в пространстве состояний

 

 

Символьный способ. Дискретную передаточную функцию получаем по формуле Q=C(z E-A)^–1B:

>>syms z;

>>А=[0 1;-6 5];В=[0; 7];С=[0 1];                   %ввод матриц

>>Q=С*inv(z*eye(2)-А))*В                            %передаточная функция

Q=7*z/(z^2-5*z+6),                                  %результат

 

Численный способ. Используем команды тулбокса CONTROL:

>> A=[0 1;-6 5];B=[0; 7];C=[0 1];

>>sd=ss(A,B,C,0,1);            % Discrete-time model Sampling time: 1,.

>> s=tf(sd)

Transfer function:

7 z ------------- z^2 - 5 z + 6

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

Моделирование дискретных систем в SIMULINK производится так же, как и непрерывных. На рисунке приведены основные блоки из библиотеки Discrete.

 

Блок задержки Unit delay позволяет строить схемы моделирования методом Кельвина. Блоки Discrete Transfer Fcn (дискретная передаточная функция) и Discrete Filter  (дискретный фильтр) определяют рациональную функцию отношения полиномов от оператора z  и обратного оператора  соответственно.

 

ЗАДАНИЕ ПО РАБОТЕ И СОДЕРЖАНИЕ ОТЧЕТА

 

Работа состоит из двух частей.

1. В первой части исследуется дискретная система, заданная разностным уравнением y_{n +2} + by_{n +1} + ay_n = u_n  с нулевыми начальными условиями y ₀=0; y ₁=0.

Требуется для своего варианта найти тремя способами ее реакцию на входной сигнал u_n =1:

· последовательно рассчитав точки y ₂, …, y ₅;

· решив разностное уравнение;

· используя z -преобразование.

Привести все числовые выкладки по решению разностного уравнения, дискретную передаточную функцию, ее разложение на простые дроби, график y_n, схему моделирования разностного уравнения в SIMULINK, построенную на элементах задержки  1/ z  и программы для вычисления в MATLAB.

2. Во второй части работы рассматривается популяция (рыб, животных), разбитых на три возрастные группы – младшего, среднего и старшего возрастов. Заданы величины р ₁, р ₂ – вероятности дожития особями каждой возрастной группы до следующего возраста, и числа а ₁, а ₂, а ₃, характеризующие среднюю плодовитость каждой возрастной группы.

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

Возрастная структура популяции описывается матричным уравнением

 

,                         (6)

 

где   x, y, z – численности трех возрастных групп.

Требуется рассчитать численности возрастных групп для n = 1, 2, 3 при начальных условиях x ₀= y ₀= z ₀=10.  Теоретически найти процентное соотношение численностей   x, y, z  для больших n – оно определяется соотношением компонент главного собственного вектора матрицы А.

Отчет должен содержать все необходимые расчеты, программы для вычисления в MATLAB и схемы моделирования в SIMULINK. Для определения собственных векторов матрицы А можно использовать символьный вариант команды e ig.

ПОРЯДОК ВЫПОЛНЕНИЯ

 

1. Получить решение разностного уравнения и построить графики

· с помощью команд ztrans, i ztrans, ezplot тулбокса SYMBOLIC;

· с помощью команд tf и step тулбокса CONTROL;

· путем моделирования в SIMULINK;

· с помощью команды rsolve пакета MAPLE.

2. Рассчитать абсолютные и относительные численности возрастных групп популяции

· с помощью степенной формулы (4); 

· с помощью команды initial.

Построить графики в обычном и логарифмическом масштабах.

КОНТРОЛЬНЫЕ ВОПРОСЫ

Даны разностные уравнения:

 

x (t +2) –5 x (t +1) + 6 x (t) = 0,	x (t +2) –3 x (t +1) + x (t) = 0,
x (t +2) – 3 x (t +1) + 2 x (t) = 0,	yk +1 + yk + 5 yk –1 + 3 yk –2 = 1,
yn +2 = yn +1 – yn,   y 1=1, y 2=0,	xk +3 + 7 xk +2 +15 xk +1 + 9 xk = 0,
xk +2–3 xk +1 + 2 xk = 0,	x (t +3)+7 x (t +2)+15 x (t +1)+9 x (t) = 0,
x (t +4) + 2 x (t +2) + x (t) = 0,	10. yn+ 3 - 5 yn+ 2 + 8 yn+ 1 - 4 yn=0;  y 0 = 0, y 1 = 2, y 2 = 10,
,	,

 

Требуется для уравнения, указанного преподавателем:

а)  Решить его с помощью характеристического уравнения;

б)  Решить его при помощи z-преобразования;

в)  Нарисовать схему моделирования и от нее перейти к матричной форме записи вида (3).

2. Решить двумя способами (по формуле (4) и путем перехода к одному уравнению) систему разностных уравнений

а)	б)	в)
Ответ

 

17.Решить систему    для трех вариантов матрицы А:

а)

Ответ для А₀:

Ответдля А₁: ^{= y ( n )}

В пределе   т.к. у матрицы А₁ главный собственный вектор [1 1 1]^T

для А₂:

В пределе при больших n .

 

ВАРИАНТЫ ЗАДАНИЙ

 

1. Варианты  коэффициентов   a и b  разностного уравнения y_{n +2} + b y_{n +1}+ a y_n =1

 

№	1	2	3	4	5	6	7	8	9	10
b	–1,3	–0,7	0,2	–2	–1,6	1,6	–1,6	1,9	–2,8	2,5
a	0,3	–0,6	–0,48	0,96	0,6	0,55	0,48	0,6	1,8	1

 

№	11	12	13	14	15	16	17	18	19	20
b	0,5	–0,72	1,28	–2	2,88	–0,5	0,72	–1,28	2	–2,88
a	0,27	0,4	0,7	1,1	1,6	0,25	0,37	0,66	1,1	1,5

 

1. Варианты  элементов  матрицы  А  системы разностных уравнений (6)  (а ₁=0)

 

№	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20
а 2	9	9 2	9	18	6	3	6	12	27 2	27 4	27 2	27	12	6	2	24	27 4	27 8	27 4	27 2
а 3	24	12	6	12	16	8	4	8	36	18	9	18	32	16	8	16	18	9	4,5	9
р 1	1 3	2 3	1 3	1 6	2 9	4 9	2 9	1 9	1 2	1	1 2	1 4	4 9	8 9	4 9	2 9	1 4	1 2	1 4	1 8
р 2	1 4	1 4	1	1	1 6	1 6	2 3	2 3	3 8	3 8	3 2	3 2	1 3	1 3	4 3	4 3	3 16	3 16	3 4	3 4

 

Литература.