Моделирование линейных непрерывных систем
При цифровом моделировании непрерывных систем необходимо обеспечить близость процессов в моделируемой непрерывной системе и в ее цифровой модели. Несовпадение этих процессов связано с двумя причинами:
1) заменой непрерывного входного процесса цифровым и 2) использованием численных методов анализа. Ошибки, связанные с заменой непрерывного процесса цифровым, были рассмотрены в предыдущей лабораторной работе. Остановимся на второй причине.
Математическая модель непрерывной системы представляет собой или нелинейное дифференциальное уравнение или совокупность соединенных между собой линейных и нелинейных блоков. В зависимости от принятой математической модели используются различные подходы к формированию цифровой модели.
Численное решение дифференциальных уравнений
Разработано большое количество методов численного решения дифференциальных уравнений. Рассмотрим, как производится численное решение на примере нелинейного дифференциального уравнения первого порядка
|
|
du / dt = f (u,x,t). (1)
Здесь x = x (t) - независимая функция (входной процесс), u = u (t) - решение уравнения (выходной процесс).
Численное решение находится для дискретных значений аргумента t, отличающихся на шаг интегрирования D t. В одношаговых разностных методах для нахождения следующего значения u к = u (t к) требуется информация только об одном предыдущем шаге. Из одношаговых методов наибольшую известность получили методы Рунге-Кутта. В основу метода Рунге-Кутта первого порядка, называемого также явным или прямым методом Эйлера, положено разложение функции u (t) в ряд Тейлора в окрестности точки A (t k-1,, u k-1):
u (t) = S0 + S1 (t - tk - 1) + S2 (t - tk - 1) 2 + …, (5.2)
где S0 = u (tk - 1) = uk - 1, Si = ( 1/ i!) du (t) / dt при t = tk - 1.
В методах Эйлера (и Рунге-Кутта тоже) ограничиваются только двумя первыми членами разложения в ряд. Запишем значение uk = u (tk), приняв в выражении (5.2) t = tk и ограничившись двумя первыми членами ряда:
uk = uk - 1 + S1 (tk - tk - 1) = uk - 1 + S1 Δ t
Учитывая, что производная du (t) / dt равна правой части дифференциального уравнения (1), имеем S1 = f (uk - 1, xk - 1, tk - 1) и окончательно получим:
uk = uk - 1 + Δ t f (uk - 1, xk - 1, tk - 1). (3)
Это выражение является приближенным решением дифференциального уравнения (1) прямым методом Эйлера. Оно рекуррентное и позволяет найти значение выходного процесса uk по значениям выходного и входного процессов в предыдущем такте.
На рис. 1 а) проиллюстрировано решение прямым методом Эйлера.
а) | б) |
Рис.1 |
Видим, что при использовании этого метода используется линейная экстраполяция и тангенс угла наклона экстраполирующей прямой равен производной функции u (t) в точке А. Экстраполированное значение uk отличается от точного на величину ошибки.
|
|
Неявный (обратный) метод Эйлера основан на разложении функции u (t) в ряд Тейлора в окрестности точки В (u k,, t k) (см. рис.1 б):
u (t) = uk + S1 (t - tk) + S2 (t - tk) 2 + …,
Приняв в этом выражении t = tk - 1 и ограничившись двумя первыми членами ряда, получим
uk - 1 = uk - Δ t f (uk, xk, tk ).
Откуда
uk = uk - 1 + Δ t f (uk, xk, tk ). (5.4)
Искомое значение процесса uk входит и в левую, и в правую части уравнения, и если не удается найти uk в явном виде, то приходится использовать приближенные методы решения этого уравнения.
Применим методы Эйлера для расчета переходной характеристики интегрирующей цепи. Передаточная функция интегрирующей цепи:
K (p) = 1/ (1 + pT).
Отсюда дифференциальное уравнение в операторной форме:
(pT + 1) y = x
и в канонической форме:
Tdy / dt + y = x.
Перепишем его в виде (1):
dy / dt = (1/ T) (x - y).
Запишем рекуррентную формулу для прямого метода Эйлера в соответствии с (5.3)
yk = yk - 1 + (Δ t / T) (xk - 1 - yk - 1), (5.5) или yk = (1 - Δ t / T) yk - 1 + Δ t / T xk - 1.
Формула для обратного метода Эйлера запишется в соответствии с (4)
yk = yk - 1 + (Δ t / T) (xk - yk).
Так как уравнение линейное, то значение yk вычисляется в явной форме:
yk = (yk - 1 + (Δ t / T) xk) / (1 + Δ t / T). (6)
Методы Эйлера обладают низкой точностью. В более точных методах используются различные способы определения угла наклона экстраполирующей прямой, чтобы она прошла ближе к точному решению. Хорошей точностью обладает метод Рунге-Кутта четвертого порядка, который обычно и используется. Программы для численного решения дифференциальных уравнений имеются практически в любом пакете прикладных программ, в том числе и в LabVIEW.
Для вычислений по формулам (5.5) и (5.6) используем структуру Formula Node. Внутри этой структуры запишем точное выражение для переходной характеристики:
z = 1 - e- i Δ t / T,
и выражения для переходной характеристики, полученные прямым методом Эйлера:
y = y 1 + (Δ t / T) (1 - y 1)
и обратным методом Эйлера:
v = (v 1 + (Δ t / T)) / (1 + Δ t / T)
при нулевых начальных условиях: y (0) = 0, v (0) = 0.
В этих выражениях использованы различные обозначения для выходных переменных и принято x = 1 (t) = 1, так как t > 0.
На рис.2 показана эта структура. В формулах Δ t обозначена как dt.
Рис.2 | Рис.3 |
Напомним, что для образования входных и выходных терминалов нужно щелкнуть ПКМ на границе структуры в предполагаемом месте терминала и в раскрывшемся меню выбрать Add Input или Add Output.
Для формирования массивов выходных переменных структура Formula Node помещается внутрь структуры For Loop, при этом задержанные на интервал дискретизации отсчеты выходных переменных y 1 и v 1 получаются с помощью регистра сдвига (рис.3).
Прямой метод Эйлера при большом интервале дискретизации может дать неустойчивое решение. Это случится, если отклонение решения от входного процесса xk - 1 - yk - 1 (см формулу (5)) даст такое значение yk. что отклонение на следующем шаге xk - yk будет той же величины, что и предыдущее, но обратным по знаку. Решение будет колебательным незатухающим.
| ||
Рис.4 |
В предыдущих лабораторных работах развертка графического индикатора Graph осуществлялась автоматически в соответствии с типом данных, подаваемых на вход графического индикатора. В этой работе мы сформируем данные так, чтобы по горизонтальной оси откладывалось время. Для этого надо сформировать кластер, куда кроме массива данных будет входить информация о времени. Используем ВП Bundle (Объединить), который находится в подпалитре Cluster (Кластер). На его входы element подаются (см. рис.4): на верхний - время начала развертки - 0; на средний - интервал дискретизации - Δt; на нижний - массив данных
|
|