Моделирование линейных непрерывных систем в среде LabVIEW

• Введение
• Моделирование линейных непрерывных систем
• Численное решение дифференциальных уравнений
• Замена непрерывной передаточной функции дискретной
• Моделирование линейных замкнутых систем
• Заключение
• Список литературы

LabVIEW (Laboratory Virtual Instrument Engineering Workbench) позволяет разрабатывать прикладное программное обеспечение для организации взаимодействия с измерительной и управляющей аппаратурой, сбора, обработки и отображения информации и результатов расчетов, а также моделирования как отдельных объектов, так и автоматизированных систем в целом. Разработчиком LabVIEW является американская компания National Instruments.

LabVIEW является открытой системой программирования и имеет встроенную поддержку всех применяемых в настоящее время программных интерфейсов, таких как Win32 DLL, COM.net, DDE, сетевых протоколов на базе IP, DataSocket и др. В состав LabVIEW входят библиотеки управления различными аппаратными средствами и интерфейсами, такими как PCI, CompactPCI/PXI, VME, VXI, GPIB (КОП), PLC, VISA, системами технического зрения и др. Программные продукты, созданные с использованием LabVIEW, могут быть дополнены фрагментами, азработанными на традиционных языках программирования, например C/С++, Pascal, Basic, FORTRAN. И наоборот можно использовать модули, разработанные в LabVIEW в проектах, создаваемых в других системах программирования. Таким образом, LabVIEW позволяет разрабатывать практически любые приложения, взаимодействующие с любыми видами аппаратных средств, поддерживаемых операционной системой компьютера.

среда программирование дифференциальное уравнение

Моделирование линейных непрерывных систем

При цифровом моделировании непрерывных систем необходимо обеспечить близость процессов в моделируемой непрерывной системе и в ее цифровой модели. Несовпадение этих процессов связано с двумя причинами:

1) заменой непрерывного входного процесса цифровым и 2) использованием численных методов анализа. Ошибки, связанные с заменой непрерывного процесса цифровым, были рассмотрены в предыдущей лабораторной работе. Остановимся на второй причине.

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

Численное решение дифференциальных уравнений

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

du/dt = f (u,x,t). (1)

Здесь x = x (t) — независимая функция (входной процесс), u = u (t) — решение уравнения (выходной процесс).

Численное решение находится для дискретных значений аргумента t, отличающихся на шаг интегрирования t. В одношаговых разностных методах для нахождения следующего значения u_к = u (t_к) требуется информация только об одном предыдущем шаге. Из одношаговых методов наибольшую известность получили методы Рунге-Кутта. В основу метода Рунге-Кутта первого порядка, называемого также явным или прямым методом Эйлера, положено разложение функции u (t) в ряд Тейлора в окрестности точки A (t_k-1,, u_k-1):

u (t) = S₀ + S₁ (t - t_{k - 1}) + S₂ (t - t_{k - 1}) ² + …, (5.2)

где S₀ = u (t_{k - 1}) = u_{k - 1,}S_i = (1/i!) du (t) /dt при t = t_{k - 1}.

В методах Эйлера (и Рунге-Кутта тоже) ограничиваются только двумя первыми членами разложения в ряд. Запишем значение u_k = u (t_k), приняв в выражении (5.2) t = t_k и ограничившись двумя первыми членами ряда:

u_k = u_{k - 1} + S₁ (t_k — t_{k - 1}) = u_{k - 1} + S₁Дt

Учитывая, что производная du (t) /dt равна правой части дифференциального уравнения (1), имеем S₁ = f (u_{k - 1}, x_{k - 1}, t_{k - 1}) и окончательно получим:

u_k = u_{k - 1} + Дt f (u_{k - 1}, x_{k - 1}, t_{k - 1}). (3)

Это выражение является приближенным решением дифференциального уравнения (1) прямым методом Эйлера. Оно рекуррентное и позволяет найти значение выходного процесса u_k по значениям выходного и входного процессов в предыдущем такте.

На рис. 1 а) проиллюстрировано решение прямым методом Эйлера.

Видим, что при использовании этого метода используется линейная экстраполяция и тангенс угла наклона экстраполирующей прямой равен производной функции u (t) в точке А. Экстраполированное значение u_k отличается от точного на величину ошибки.

Неявный (обратный) метод Эйлера основан на разложении функции u (t) в ряд Тейлора в окрестности точки В (u_k,, t_k) (см. рис. 1 б):

u (t) = u_k + S₁ (t - t_k) + S₂ (t - t_k) ² + …,

Приняв в этом выражении t = t_{k - 1} и ограничившись двумя первыми членами ряда, получим

u_{k - 1} = u_k — Дt f (u_k, x_k, t_k).

Откуда

u_k = u_{k - 1} + Дt f (u_k, x_k, t_k). (5.4)

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

Применим методы Эйлера для расчета переходной характеристики интегрирующей цепи. Передаточная функция интегрирующей цепи:

K (p) = 1/ (1 + pT).

Отсюда дифференциальное уравнение в операторной форме:

(pT + 1) y = x

и в канонической форме:

Tdy/dt + y = x.

Перепишем его в виде (1):

dy/dt = (1/T) (x - y).

Запишем рекуррентную формулу для прямого метода Эйлера в соответствии с (5.3)

y_k = y_{k - 1} + (Дt/T) (x_{k - 1} — y_{k - 1}), (5.5) или y_k = (1 — Дt/T) y_{k - 1} + Дt/T x_{k - 1}.

Формула для обратного метода Эйлера запишется в соответствии с (4)

y_k = y_{k - 1} + (Дt/T) (x_k — y_k).

Так как уравнение линейное, то значение y_k вычисляется в явной форме:

y_k = (y_{k - 1} + (Дt/T) x_k) / (1 + Дt/T). (6)

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

Для вычислений по формулам (5.5) и (5.6) используем структуру Formula Node. Внутри этой структуры запишем точное выражение для переходной характеристики:

z= 1 — e^-iДt/T,

и выражения для переходной характеристики, полученные прямым методом Эйлера:

y = y₁ + (Дt/T) (1 — y₁)

и обратным методом Эйлера:

v = (v₁ + (Дt/T)) / (1 + Дt/T)

при нулевых начальных условиях: y (0) = 0, v (0) = 0.

В этих выражениях использованы различные обозначения для выходных переменных и принято x = 1 (t) = 1, так как t > 0.

На рис. 2 показана эта структура. В формулах Дt обозначена как dt.

Напомним, что для образования входных и выходных терминалов нужно щелкнуть ПКМ на границе структуры в предполагаемом месте терминала и в раскрывшемся меню выбрать Add Input или Add Output.

Для формирования массивов выходных переменных структура Formula Node помещается внутрь структуры For Loop, при этом задержанные на интервал дискретизации отсчеты выходных переменных y₁ и v₁ получаются с помощью регистра сдвига (рис.3).

Прямой метод Эйлера при большом интервале дискретизации может дать неустойчивое решение. Это случится, если отклонение решения от входного процесса x_{k - 1} — y_{k - 1} (см формулу (5)) даст такое значение y_k. что отклонение на следующем шаге x_k — y_k будет той же величины, что и предыдущее, но обратным по знаку. Решение будет колебательным незатухающим.

В предыдущих лабораторных работах развертка графического индикатора Graph осуществлялась автоматически в соответствии с типом данных, подаваемых на вход графического индикатора. В этой работе мы сформируем данные так, чтобы по горизонтальной оси откладывалось время. Для этого надо сформировать кластер, куда кроме массива данных будет входить информация о времени. Используем ВП Bundle (Объединить), который находится в подпалитре Cluster (Кластер). На его входы element подаются (см. рис.4): на верхний — время начала развертки — 0; на средний — интервал дискретизации — Дt; на нижний — массив данных

Замена непрерывной передаточной функции дискретной

Если математическая модель системы представляется в виде соединения линейных и нелинейных блоков, то для описания линейных блоков чаще всего используется передаточная функция K (p). В этом случае цифровую модель непрерывного линейного блока можно получить, заменив непрерывную передаточную функцию K (p) дискретной K (z).

Для этого можно использовать связь между непрерывными и дискретными изображениями, устанавливаемую дискретным преобразованием Лапласа (Z-преобразованием). В таблице 1 приведена эта связь для передаточных функций, используемых в данной лабораторной работе.

Таблица 1

Заметим, что здесь комплексная переменная z определяется как z = e^pДt и является оператором опережения на интервал дискретизации. Соответственно z^-1 — это оператор задержки на интервал дискретизации.

Другой путь предусматривает непосредственный переход от комплексной переменной p к комплексной переменной z заменой операции аналогового интегрирования 1/p операцией дискретного интегрирования. При дискретном описании аналогового интегрирования можно оперировать только с значениями входного и выходного процессов в моменты дискретизации. На рис. 5 показано, как это можно сделать, используя численное интегрирование по методу прямоугольников и по методу трапеций.

Значение выходного процесса y_k интегратора в момент времени t = kДt отличается от предыдущего значения y_k-1 на величину площади S под кривой x (t) (заштрихованная фигура на рис. 5 а).

По методу прямоугольников площадь можно определить по разному в зависимости от того, какую величину принять за высоту прямоугольника: x_k-1 или x_k (рис. 5.5 б и рис. 5.5 в). На рис. 5.5 г) показано, как вычисляется эта площадь по методу трапеций. Рекуррентные формулы для интегрирования приведены под рисунками.

По этим формулам можно записать дискретные передаточные функции. Поясним это на примере интегрирования по методу трапеций:

y_k = y_k-1 + Дt (x_k + x_k-1) /2.

Перенесем y_k-1 в левую часть и возьмем от полученного выражения Z-преобразование. Учитывая, что запаздывание на интервал дискретизации в области оригиналов соответствует умножению на z^-1 в области изображений, получим:

Y (z) — z^-1Y (z) = (Дt/2) (X (z) + z^-1X (z)).

Дискретная передаточная функция — это отношение Z-изображений выходной и входной переменных, поэтому

K (z) = Y (z) /X (z) = (Дt/2) (1 + z^-1) / (1 — z^-1) = (Дt/2) (z + 1) / (z — 1).

В таблице 2 приведены выражения дискретных передаточных функций для различных методов численного интегрирования для одного и двух интеграторов.

Таблица 2

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

В таблице 1 была приведена дискретная передаточная функция интегрирующей цепи (для которой К (р) = 1/ (1 + рТ)), полученной применением Z-преобразования. Найдем другие варианты дискретной передаточной функции интегрирующей цепи, отличающиеся методами численного интегрирования.

При использовании метода прямоугольников (1) в передаточную функцию K (p) = 1/ (1 + pT) вместо р нужно подставить (z — 1) /Дt. Тогда получим

K (z) = 1/ (1 + (z — 1) T/Дt) =.

Аналогично можно получить дискретные передаточные функции и для других методов численного интегрирования. Они представлены в таблице 3 Принято обозначение Дt/T = б Таблица 3

Этим передаточным функциям соответствуют следующие рекуррентные формулы.

Для Z-преобразования

y_k = e^-бy_{k - 1} + бx_k. (7)

Для численного интегрирования по методу прямоугольников (1)

y_k = (1 — б) y_{k - 1} + бx_{k - 1}.

Полученная формула совпадает с формулой для прямого метода Эйлера Для численного интегрирования по методу прямоугольников (2)

y_k = (1/ (1 + б)) y_{k - 1} + (б/ (1 + б)) x_k. (8)

и по методу трапеций

y_k = ((2 — б) / (2 + б)) y_{k - 1} + (б/ (2 + б)) (x_k + x_{k - 1}). (9)

В лабораторной работе производится оценка ошибок цифрового моделирования для каждого из этих методов.

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

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

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

Проводя аналогию с теоремой Котельникова, можно потребовать, чтобы Дf_0,1Дt = 5 — 10, где Дf_0,1 — полоса пропускания замкнутой системы по уровню 0,1.

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

К_р (р) =. (10)

Такая модель часто используется при анализе ошибок в следящей системе.

Запишем дискретную передаточную функцию разомкнутой системы, заменяя интеграторы по методу прямоугольников (2). Для этого преобразуем передаточную функцию разомкнутой системы (5.10), поделив числитель и знаменатель на р²:

К_р (р) =.

Используя соотношения, приведенные в таблице (5.2), получим:

К_р (z) = =, (11)

где b₁ = K (Дt) ²/ (T + Дt), a₁ = - (2T + Дt) / (T + Дt), a₂ = T/ (T + Дt).

Для моделирования устройства с передаточной функцией (11) используется БИХ-фильтр, коэффициенты числителя которого (Forward Coefficients) представляются массивом из двух элементов (0,b₁), а коэффициенты знаменателя — массивом из трех элементов (1, a₁,a₂).

Специфика использования БИХ-фильтра заключается в том, что неизвестен целиком входной массив Х, а известен только текущий элемент, а следующий элемент рассчитывается с учетом значения текущего элемента выходного массива фильтра. В LabVIEW существует такой фильтр — IIR Filter PtByPt (IIR Filter Point By Point — БИХ-фильтр точка за точкой).

Вычисления БИХ-фильтром IIR Filter PtByPt производятся в цикле For

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

Для точного расчета переходной характеристики воспользуемся ВП ODE Linear nth Order Numeric — «Решение линейного обыкновенного дифференциального уравнения n-го порядка в численном виде» (рис.7).

ВП находит решение в виде суммы экспонент и вычисляет его для заданных точек. Поэтому решение точное.

Вход, А представляет собой массив коэффициентов дифференциального уравнения в порядке увеличения степени производной. Коэффициент при производной самой высокой степени считается равным 1 и не требует ввода.

На вход Х0 подается массив начальных условий — начальные значения решения и его n — 1 — й производных.

Вход «число точек» задает число равноудаленных по времени точек между начальным и конечным временем

Выход Х содержит массив значений решения в равномерно расположенных по оси времени точках. Значение времени в этих точках выводится в массиве Times.

Дифференциальное уравнение замкнутой системы запишем по передаточной функции замкнутой системы:

К_з (р) = = (12)

Дифференциальное уравнение замкнутой системы

Td²y/dt² + dy/dt + Ky = Kx

Запишем однородное дифференциальное уравнение, учитывая, что коэффициент при высшей производной должен быть равен 1

d²y/dt² + (1/Т) dy/dt + (K/Т) y = 0

Для компьютерного решения этого уравнения нужно задать массив, А = (К/Т, 1/Т). Чтобы получить переходную характеристику, нужно задать массив Х = (- 1, 0) и к решению прибавить 1.

Полностью блок-схема программы моделирования замкнутой системы приведена на рис. 8.

Рис.8

Заключение

В основе технологии использования LabVIEW лежит комбинированное моделирование систем на ЭВМ, включающее аналитическое, имитационное и натурное.

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

Натурным моделированием называют проведение исследования на реальном объекте с возможностью вмешательства человека в процесс проведения эксперимента и последующей обработки результатов эксперимента на вычислительной технике.

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

1. Н. А. Виноградова, Я. И. Листратов, Е. В. Свиридов. «Разработка прикладного программного обеспечения в среде LabVIEW». Учебное пособие — М.: Издательство МЭИ, 2005.

2. http://www.automationlabs.ru/

3. http://digital. ni.com/

4. http://www.labview.ru/

5. http://ru. wikipedia.org/