Среди способов моделирования непрерывных систем управления электроприводом можно выделить два, основанных на использовании математических моделей систем в виде моделей состояния и структурных моделей, каждый из которых имеет свои определенные преимущества при решении конкретных задач моделирования АСУ ЭП. Наиболее удобно использовать модель состояния при моделировании и синтезе многомерных линейных систем управления ЭП методами пространств состояний. При моделировании нелинейных систем ЭП, а также некоторых специфических элементов современных систем ЭП, например тиристорных преобразователей и микропроцессоров, более эффективным является использование структурных моделей. Особенно удобно их применять при анализе в связи с выраженной структурой реальных систем электропривода. Однако эффективность использования структурных (топологических) методов существенно снижается по мере усложнения систем управления ЭП. Поэтому выбор способа моделирования обуславливается целесообразностью его применения в конкретном случае.

Цифровое моделирование непрерывных систем управления основывается на описании системы обыкновенными дифференциальными уравнениями в форме Коши, где в общем случае для многомерного элемента каждая входная переменная связана с каждой выходной переменной. Если взаимосвязи по всем каналам линейны или линеаризованы, то в общем случае многомерный элемент можно описать системой неоднородных дифференциальных уравнений. Систему можно записать более компактно в виде одного векторного дифференциального уравнения. Векторное дифференциальное уравнение в форме Коши, отражающее динамические свойства многомерного линейного объекта, является уравнением состояния и используется в качестве математической модели при моделировании методами пространств состояний. Полная математическая модель линейного многомерного объекта, кроме уравнений состояния, содержит еще уравнение выхода, связывающее переменные состояния и управляющие воздействия с выходными переменными.

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

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

Матричные методы расчета переходного процесса в линейных системах основаны на расчете переходной (экспоненциальной) матрицы состояния, что связано с необходимостью выполнения сложных и громоздких расчетов, и особенно затруднены при отсутствии специализированных прикладных пакетов программ (наиболее известным пакетом символьной математики, ориентированным на работу с векторами и матрицами, следует признать MatLab). Способы вычисления переходной матрицы состояния можно классифицировать следующим образом: прямые, основанные на методике Планта, аппроксимации Паде, теоремы Кели-Гамильтона. Все перечисленные методы вычисления переходной матрицы состояния используют рекуррентный алгоритм ее расчета. Переходная матрица состояния представляется разложением в матричный ряд. Для обеспечения работоспособности алгоритма вычисления переходной матрицы необходимо установить максимальное число членов ряда, при превышении которого вычисления прекращаются. Следует отметить, что при числе членов ряда к =2 точность вычисления переходной матрицы состояния соответствует точности метода Эйлера, при к =3 - точности усовершенствованного метода Эйлера, при к =5 - точности метода Рунге-Кутта. Очевидно, что затраты на вычисления значительно выше по сравнению с методами численного интегрирования. Кроме выполнения расчетов для переходной матрицы состояния, необходимо выполнить вычисление входной матрицы, при котором используются в основном два метода: аналитический, когда заранее известно, что переходный процесс имеет устойчивый характер; приближенный, когда характер переходного процесса заранее не определен. Использование и того и другого метода связано с громоздкими матричными операциями. Но следует отметить, что матричный метод имеет свои преимущества перед остальными методами при моделировании многомерных систем управления, имеющих несколько входов и выходов.

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

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

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

Результаты изучения переходных процессов в физических системах на основе методов, в которых непрерывные сигналы при расчете заменяются временными последовательностями чисел, показывают, что такой подход дает значительную экономию вычислительных затрат. Соотношения между временными последовательностями действительных чисел (решетчатыми функциями) описываются удобными рекуррентными разностными уравнениями, коэффициенты которых зависят от параметров физических систем. Некоторые рекуррентные методы, в частности метод Тастина, позволяют получать эффективные алгоритмы цифрового моделирования дискретных систем. Сущность известных в настоящее время рекуррентных разностных методов и состоит в замене процессов, происходящих в непрерывных системах, процессами в эквивалентных дискретных системах. Математическим аппаратом при этом служит метод z -преобразований. Рассмотренные методы Тастина, Боксера-Талера построения цифровых моделирующих алгоритмов систем управления, заданных в виде структурных схем, имеют очень мало ограничений или вообще не имеют. Они являются универсальными в смысле использования при входных сигналах аналитической или произвольной формы. Порядок рекуррентных уравнений совпадает с порядком линейной части моделируемой системы независимо от используемого метода. Не требуется дополнительных усилий при проведении подготовительной работы. Однако точность этих методов принципиально не так высока, как методов использующих информацию о всей непрерывной системе в целом (методы инвариантных импульсных функций, Цыпкина-Гольденберга, Рагаццини-Бергена).

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

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

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

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

Одним из основных параметров при имитационном моделиро­вании является модельное время, которое отображает время функ­ционирования реальной системы. В зависимости от способа продвижения модельного времени методы моделирования подраз­деляются на методы с приращением временного интервала и методы с продвижением времени до особых состояний. В первом случае модельное время продвигается на некоторую величину Dt . Определяются изменения состояний элементов и выходных воздей­ствий системы, которые произошли за это время. После этого модельное время снова продвигается на величину Dt , и процедура повторяется. Так продолжается до конца периода моделирова­ния Т m . Шаг приращения времени Dt зачастую выбирается по­стоянным, но в общем случае он может быть и переменным. Этот метод называют «принципом Dt ».

Во втором случае в текущий момент модельного времени t сначала анализируются те будущие особые состояния – поступле­ние дискретного входного воздействия (заявки), завершение обслуживания и т.п., для которых определены моменты их на­ступления t i > t. Выбирается наиболее раннее особое состояние, и модельное время продвигается до момента наступления этого состояния. Считается, что состояние системы не изменяется между двумя соседними особыми состояниями. Затем анализируется реакция системы на выбранное особое состояние. В частности, в ходе анализа определяется момент наступления нового особого состояния. Затем анализируются будущие особые состояния, и модельное время продвигается до ближайшего. Процедура повто­ряется до завершения периода моделирования T m . Данный метод называют «принципом особых состояний», или «принципом dz ». Благодаря его применению экономится машинное время моделирования. Однако он используется только тогда, когда имеется возможность определения моментов наступления будущих очередных особых состояний.

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

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

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

Общие сведения

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

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

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

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

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

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

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

u (t ) = S 0 + S 1 (t – t k - 1) + S 2 (t – t k - 1) 2 + …, (5.2)

где S 0 = 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 1 (t k – t k - 1) = u k - 1 + S 1 Δt

Учитывая, что производная du (t )/dt равна правой части дифференциального уравнения (5.1), имеем S 1 = 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). (5.3)

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

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

а )	б )
Рис. 5.1

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

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

u (t ) = u k + S 1 (t – t k ) + S 2 (t – t k ) 2 + …,

Приняв в этом выражении 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 .

Перепишем его в виде (5.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 .

Формула для обратного метода Эйлера запишется в соответствии с (5,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 ). (5. 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.

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

Рис. 5.2	Рис. 5.3

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

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

Прямой метод Эйлера при большом интервале дискретизации может дать неустойчивое решение. Это случится, если отклонение решения от входного процесса x k – 1 – y k – 1 (см формулу (5.5)) даст такое значение y k . что отклонение на следующем шаге x k – y k будет той же величины, что и предыдущее, но обратным по знаку. Решение будет колебательным незатухающим.

К графическому индикатору

Рис. 5.4

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

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

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

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

Таблица 5.1

K (p )	1 p	p 2	(1 + pT )
K (z )	Δt z (z – 1)	(Δt ) 2 z (z – 1) 2	(Δt /T )z (z – e - Δt /T )

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

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

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

y k = y k- 1 + S	y k = y k- 1 + Δt x k- 1	y k = y k- 1 + Δt x k	y k = y k - 1 + +Δt (x k + x k - 1)/2
а )	б )	в )	г )
Рис. 5.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 - 1 Y (z ) = (Δt /2)(X (z ) + z - 1 X (z )).

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

K (z ) = Y (z )/X (z ) = (Δt /2)(1 + z -1)/(1 – z -1) = (Δt /2)(z + 1)/(z – 1).

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

Таблица 5.2

K (p )	K (z )
Метод прямоугольников (1)	Метод прямоугольников (2)	Метод трапеций
p	Δt z – 1	Δt z z – 1	Δt (z + 1) 2 (z – 1)
p 2	(Δt ) 2 (z +1) 2(z – 1) 2	(Δt ) 2 z (z – 1) 2	(Δt ) 2 (z 2 + 4z + 1) 6(z – 1) 2

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

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

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

	Δt /T z – (1 – Δt /T )

K (z ) = 1/(1 + (z – 1)T /Δt ) = .

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

Таблица 5.3

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

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

y k = e - α y k – 1 + αx k . (5.7)

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

y k = (1 – α)y k – 1 + αx k - 1 .

Полученная формула совпадает с формулой для прямого метода Эйлера

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

y k = (1/(1 + α))y k – 1 + (α/(1 + α))x k . (5.8)

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

y k = ((2 – α)/(2 + α))y k – 1 + (α/(2 + α))(x k + x k – 1). (5.9)

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

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

Нужно быть очень внимательным при выборе интервала дискретизации, когда моделируются замкнутые системы. В этих системах текущее значение входного процесса сравнивается со значением выходного процесса, рассчитанного по предыдущим значениям входного процесса. Это экстраполированное значение не должно значительно отличаться от входного процесса. В противном случае возникают большие ошибки моделирования, а при большом интервале дискретизации процесс может стать неустойчивым. Выбор интервала дискретизации нужно связывать с полосой пропускания замкнутой системы. Проводя аналогию с теоремой Котельникова, можно потребовать, чтобы Δf 0,1 Δt = 5 – 10, где Δf 0,1 – полоса пропускания замкнутой системы по уровню 0,1.

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

К р (р ) = .

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

	К (Δt ) 2 z /(z – 1) 2 Т + Δt z /(z – 1)		b 1 z z 2 + a 1 z + a 2

К р (z ) = = , (5.11)

где b 1 = K (Δt ) 2 /(T + Δt ), a 1 = - (2T + Δt )/(T + Δt ), a 2 = T /(T + Δt ).

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

Рис. 5.6

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

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

Рис. 5.7

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

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

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

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

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

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

	К р (р ) 1 + К р (р )		К Тр 2 + р + К

К з (р ) = = , (5.12)

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

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

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

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

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

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

k K T k dt

Рис. 5.8

Контрольные вопросы

1. В каком виде записывается нелинейное дифференциальное уравнение первого порядка для численного решения?

2. Как записывается разложение функции в ряд Тейлора?

3. Поясните графически решение дифференциального уравнения прямым методом Эйлера.

4. Как записывается рекуррентная формула для численного решения нелинейного дифференциального уравнения первого порядка прямым методом Эйлера?

5. Поясните графически решение дифференциального уравнения обратным методом Эйлера.

6. Как записывается рекуррентная формула для численного решения нелинейного дифференциального уравнения первого порядка обратным методом Эйлера?

7. Запишите дифференциальное уравнение интегрирующей цепи в форме, удобной для решения методом Эйлера.

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

9. Запишите рекуррентную формулу для численного решения дифференциального уравнения интегрирующей цепи обратным методом Эйлера.

10. Какие два пути используются при замене непрерывной передаточной функции дискретной передаточной функцией?

11. Запишите рекуррентную формулу для численного интегрирования по методу прямоугольников (1).

12. Запишите рекуррентную формулу для численного интегрирования по методу прямоугольников (2).

13. Запишите рекуррентную формулу для численного интегрирования по методу трапеций.

14. Выведите дискретную передаточную функцию интегратора по методу прямоугольников (1).

15. Выведите дискретную передаточную функцию интегратора по методу прямоугольников (2).

16. Выведите дискретную передаточную функцию интегратора по методу трапеций.

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

18. Какая структура используется для вычислений по рекуррентным формулам?

19. Откуда при моделировании берется значение y k – 1 , необходимое для расчета y k ?

20. Как образовать терминалы для ввода и вывода переменных в структуре Formula Node?

21. Почему при моделировании замкнутой системы используется ВП IIR Filter PtByPt, а не ВП IIR Filter?

22. Почему при соединении выхода БИХ-фильтра с его входом в цепи обратной связи автоматически появляется регистр сдвига?

23. Для чего используется ВП ODE Linear nth Order Numeric?

24. Для чего используется ВП Bundle (Объединить)?

Программа работы

1.Вызвать пакет LabVIEW. Открыть New VI.

2. Численное решение дифференциального уравнения интегрирующей цепи методами Эйлера.

2.1. Сформировать лицевую панель создаваемого ВП. На лицевую панель вывести цифровые элементы управления “Постоянная времени Т”, “Интервал дискретизации dt”, “Множитель развертки k”. Поместить на лицевую панель графический индикатор Graph. Разместить под графическим индикатором палитру элементов управления графиком (Graph Palette) и панель редактирования курсора (вывести два курсора разного цвета).

2.2. В окно BD поместить структуру Formula Node. Поместить внутри нее программу расчета переходной характеристики интегрирующей цепи по точной формуле и по формулам для прямого и обратного методов Эйлера (рис. 5.2). На границе структуры поместить терминалы для всех входных и выходных величин. Размер структуры взять таким, чтобы в ней можно было разместить еще четыре строчки.

2.3. Поместить структуру Formula Node внутрь структуры For Loop. Задать число итераций равным (T*k/dt). На границах структуры For Loop поместить терминалы регистров сдвига для переменных y и v, чтобы получить задержанные значения y1 и v1. Установить нулевые начальные условия (рис. 5.3).

2.4. Подсоединить к терминалам dt и Т структуры Formula Node выходы соответствующих элементов управления.

2.5. Одномерные массивы выходных величин z, y и v объединить в многомерный массив, используя функцию Build Array.

2.6. Сформировать кластер, используя функцию Bundle (рис.5.4). К верхнему входу element подсоединить 0, к среднему – dt, к нижнему массив с выхода ВП Build Array. Убедиться в отсутствии ошибок в собранной блок-диаграмме.

2.7. Задать постоянную времени в соответствии с выполняемым вариантом.

Вариант
Т , с	0.1	0.4	0.8	0,2	0,5	0,9	0,3	0,6	1,0	0,7

Интервал дискретизации dt взять в 20 раз меньше постоянной времени Множитель развертки – 5 – 10. Запустить моделирование. Убедиться в близости всех решений.

2.8 Исследовать влияние интервала дискретизации на ошибку решения дифференциального уравнения методами Эйлера. Вывести вертикальную линию курсора на время, равное 2Т. Закрепить один курсор на точном решении (белый график): щелчок ЛКМ на кнопке в панели редактирования курсоров. В раскрывшемся меню выбрать Plot 0. Второй курсор поочередно закреплять на Plot 1 и Plot 2. Проделать замеры ошибок для прямого и обратного методов Эйлера для интервалов дискретизации равных 0,05Т; 0,1Т; 0,2Т; 0,5Т; Т; 2Т. При каком интервале дискретизации решение прямым методом Эйлера становится неустойчивым? Объясните, почему.

3. Моделирование интегрирующей цепи заменой непрерывной передаточной функции дискретной передаточной функцией

3.1. Поместить Разработанную блок-диаграмму в структуру Case Structure так, чтобы вся блок-диаграмма кроме элементов управления и графического индикатора с формирователем кластера оказалась внутри структуры. Напомним работу с этой структурой. Переключение вариантов произведем с помощью строковой переменной. Активизировать FP. Поместить на лицевую панель переключатель вариантов: Controls → All Controls → Ring & Enum → Enum. Озаглавить его “Варианты”. Активизировать BD. Подсоединить выход узла “Варианты” к терминалу селектора вариантов . В переключателе вариантов (на верхней границе структуры) вместо True и False появятся цифры 1 и 0. Заметим, что собранный в п.2 генератор соответствует цифре 1. Назвать варианты: щелкнуть ПКМ на узле Enum; в появившемся меню выбрать опцию Properties. В раскрывшемся окне Enum Properties открыть закладку Edit Items. В области Items набрать названия вариантов: “Передаточные функции” (ему присваивается символ 0) и “Дифференциальное уравнение” (символ 1). Подтвердить установку – ОК.

3.2. Скопировать всю схему варианта “Дифференциальное уравнение” и вставить ее в вариант “Передаточные функции”. Восстановить все нарушенные соединения.

3.3. Моделирование будем проводить по рекуррентным формулам (5.7) – (5.9). Перепишем их, обозначив переменные другими буквами.

y k = e - α y k – 1 + αx k ,

v k = (1/(1 + α))v k – 1 + (α/(1 + α))x k ,

w k = ((2 – α)/(2 + α))w k – 1 + (α/(2 + α))(x k + x k – 1)

Напомним, что α = dt/T.

Внести изменения в записи внутри структуры Formula Node в соответствии с рис. 5.9.

Рис. 5.9

Поместить терминалы на границе структуры для вновь введенных переменных. Ввести для новых переменных регистры сдвига в структуре For Loop.

Четыре выходных массива z, y, v и w объединить функцией Build Array. Убедиться в отсутствии ошибок.

3.4. Замерить ошибки для этих трех методов, используя курсоры для того же значения постоянной времени Т, что и в п.2 для значений интервала дискретизации dt, равного 0,1Т и 0,5Т для того же момента времени 2Т.

3.5. Сравнить исследованные в п.2 и п.3 методы по точности, приведя необходимые графики.

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

4.1. Добавить в переключатель вариантов “Варианты” к имеющимся двум еще один, назвав его “Замкнутая система”. Для этого щелкнуть ПКМ по узлу Варианты (Enum), в раскрывшемся меню выбрать Properties, в окне Enum Properties открыть закладку Edit Items. В области Items набрать указанное название варианта. Ему автоматически присвоится цифра 2. Щелкнуть ПКМ по границе структуры Case. В раскрывшемся меню выбрать Add Case for Every Value (Добавить варианты для каждого значения).

4.2. На лицевой панели поместить цифровой элемент управления, назвав его “Коэффициент передачи”

4.3. Собрать блок-схему программы по рис. 5.8. Количество итераций N такое же, как и в предыдущих вариантах.

4.6. Убедиться в отсутствии ошибок. Задать значения К и Т в соответствии с выполняемым вариантом

Вариант
К
Т, с	0,05	0,033	0,02	0,012	0,008	0,04	0,025	0,017	0,01	0,007

4.7. Исследовать влияние интервала дискретизации на ошибку моделирования приняв за показатель точности величину первого выброса в переходной характеристике h m , а отношение Dh m /h m – за относительную ошибку моделирования. Здесь Dh m – разность между h m при выбранном шаге моделирования и точном значении h m . За точное значение h m принять его значение при Δt = 0,0001

Результаты измерений свести в таблицу:

Лабораторная работа №6

Методы моделирования систем

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

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

При этом возможно изменение не только модели, но и метода моделирования, что часто является средством развития представления ЛПР о моделируемой ситуации. Иными словами, перевод вербального описания в формальное, осмысление, интерпретация модели и получаемых результатов становятся неотъемлемой частью практически каждого этапа моделирования сложной развивающейся системы.

Часто для того чтобы точнее охарактеризовать такой подход к моделированию процессов принятия решений, говорят о создании «механизма» моделирования, «механизма» принятия решений (например, «хозяйственный механизм», «механизм проектирования и развития предприятия» и т.п.).

Возникающие вопросы – как формировать такие развивающиеся модели или «механизмы»? как доказывать адекватность моделей? – и являются основным предметом системного анализа.

Для решения проблемы перевода вербального описания в формальное в различных областях деятельности стали развиваться специальные приёмы и методы. Так, возникли методы типа «мозговой атаки», «сценариев», экспертных оценок, «дерева целей» и т.п.

В свою очередь, развитие математики шло по пути расширения средств постановки и решения трудноформализуемых задач. Наряду с детерминированными, аналитическими методами классической математики возникла теория вероятностей и математическая статистика (как средство доказательства адекватности модели на основе представительной выборки и понятия вероятности правомерности использования модели и результатов моделирования). Для задач с большей степенью неопределённости инженеры стали привлекать теорию множеств, математическую логику, математическую лингвистику, теорию графов, что во многом стимулировало развитие этих направлений. Иными словами, математика стала постепенно накапливать средства работы с неопределённостью, со смыслом, который классическая математика исключала из объектов своего рассмотрения.

Таким образом, между неформальным, образным мышлением человека и формальными моделями классической математики сложился как бы «спектр» методов, которые помогают получать и уточнять (формализовать) вербальное описание проблемной ситуации, с одной стороны, и интерпретировать формальные модели, связывать их с реальной действительностью, с другой. Этот спектр условно представлен на рис. 2.1, а.

Развитие методов моделирования, разумеется, шло не так последовательно, как показано на рис. 2.1, а. Методы возникали и развивались параллельно. Существуют различные модификации сходных методов. Их по-разному объединяли в группы, т.е. исследователи предлагали разные классификации (в основном – для формальных методов, что более подробно будет рассмотрено в следующем параграфе). Постоянно возникают новые методы моделирования как бы на «пересечении» уже сложившихся групп. Однако основную идею – существование «спектра» методов между вербальным и формальным представлением проблемной ситуации – этот рисунок иллюстрирует.

Первоначально исследователи, развивающие теорию систем, предлагали классификации систем и старались поставить им в соответствие определённые методы моделирования, позволяющие наилучшим образом отразить особенности того или иного класса. Такой подход к выбору методов моделирования подобен подходу прикладной математики. Однако в отличие от последней, в основу которой положены классы прикладных задач, системный анализ может один и тот же объект или одну и ту же проблемную ситуацию (в зависимости от степени неопределённости и по мере познания) отображать разными классами систем и соответственно различными моделями, как бы организовывая таким образом процесс постепенной формализации задачи, т.е. «выращивание» её формальной модели. Подход помогает понять, что неверно выбранный метод моделирования может привести к неверным результатам, к невозможности доказательства адекватности модели, к увеличению числа итераций и затягиванию решения проблемы.

Лабораторная работа №2

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

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

При замене непрерывного процесса цифровым возникают ошибки из-за квантования по уровню (шумы квантования) и дискретизации по времени (ошибки восстановления непрерывного процесса по его дискретным отсчетам). Шумы квантования считаются случайным процессом с дисперсией h 2 /12, где h – величина шага квантования. При 16-разрядном двоичном представлении числа шаг квантования равен примерно одной 65-тысячной этого числа. Поэтому шумами квантования можно пренебречь. Восстановить непрерывный процесс по его дискретным отсчетам можно без ошибки согласно теореме Котельникова, если спектр этого процесса ограничен частотой f гр (S (f ) = 0 при f > f гр) и частота дискретизации f дискр ³ 2f гр. При функциональном моделировании систем частоту дискретизации обычно берут много больше f гр: f дискр = (10 – 20)f гр, где f гр = f в – верхняя частота спектра, т. е. частота, на которой спектр процесса спадает до достаточно малой величины.

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

Разработано большое количество методов численного решения дифференциальных уравнений. Рассмотрим, как производится численное решение на примере нелинейного дифференциального уравнения первого порядка вида du /dt = f (u , t ). Решение находится для дискретных значений аргумента, отличающихся на шаг интегрирования Dt . В одношаговых разностных методах для нахождения следующего значения u к = u (t к) требуется информация только об одном предыдущем шаге. Из одношаговых методов наибольшую известность получили методы Рунге-Кутта. В основу метода Рунге-Кутта первого порядка, называемого также явным методом Эйлера, положено разложение функции u (t ) в ряд Тейлора в окрестности точки (t k -1, u k -1), ограниченное двумя первыми членами ряда: u k = u k – 1 + Dt *u " k -1 , где u " k -1 = du (t )/dt при t = t r -1 . Так как du (t )/dt = f (u , t ), то u k =u k – 1 +Dt *f (u k – 1 ,t k –1). Видим, что при использовании этого метода считается, что в течение времени Dt функция u (t ) изменяется линейно и тангенс угла наклона прямой равен u " k -1 . Это, как показано ниже на рисунке, приводит к ошибке (сплошные линии).

Неявный (обратный) метод Эйлера основан на разложении функции u (t ) в ряд Тейлора в окрестности точки (t k , u k). Расчет ведется по выражению u k = u k – 1 + Dt *u " k = u k – 1 + Dt *f (u k ,t k). Решение находится тоже с ошибкой, хотя и другого знака (пунктирная линия на
рисунке). Ошибка увеличивается с увеличением шага Dt = t k - t k -1 . Для уменьшения этой ошибки при неизменном Dt используют методы Рунге-Кутта более высокого порядка. При большом шаге вычислительный процесс может стать неустойчивым.

VisSim является пакетом визуального моделирования, поэтому задать модель аналитически в виде нелинейного дифференциального уравнения нельзя. Нужно по дифференциальному уравнению составить функциональную схему, или графическую модель. Для дифференциального уравнения du /dt = f (u , t ) она составляется следующим образом. Сначала формируется f (u , t ) нелинейным блоком, на входы которого подаются t и u . Так как f (u , t ) является производной процесса u (t ), то сам процесс u (t ) получается интегрированием выходного процесса нелинейного блока. Получившаяся графическая модель показана на рисунке ниже.

Для линейных систем возможен другой путь расчета выходного процесса – с использованием интеграла свертки u вых (t ) = ∫ u вх (t - τ)g (τ)d τ. Для расчета u k используются численные методы вычисления интеграла с верхним пределом t = k Δt . В пакете VisSim именно он используется для расчета процесса на выходе линейных звеньев, задаваемых передаточными функциями. При этом рассчитанное значение выходного процесса в любой момент времени при скачкообразном входном процессе будет находиться на одной и той же линии переходной характеристики. От шага модельного времени зависит только количество рассчитанных точек переходной характеристики. Если входной процесс произволен, то шаг модельного времени в этом случае определяется из соображений допустимой ошибки при дискретизации процесса.

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

Выполнение работы

1. Исследование ошибок численного решения дифференциального уравнения.

1.1 Составить графическую модель для дифференциального уравнения согласно заданию:

Уравнение	du /dt = at 2 + bu			du /dt = at (1 + bu )			du /dt = at 2 + but

1.2. Собрать на рабочем столе VisSim составленную модель. Запомните, что различные методы численного интегрирования реализуются только для моделей, в которые входят интеграторы, обозначаемые как 1/S (Blocks → Integration → integrator). Интегратор можно реализовать и по-другому – как линейное устройство, задаваемое передаточной функцией (Blocks → Linear System → transferFunction), если установить в числителе (Numenator) – 1, а в знаменателе (Denominator) – 1_0. Но для такого интегратора, охваченного обратной связью, реализуется только метод Эйлера. Для удобства изображения линии обратной связи использовать указатель (Blocks → Annotation → wirePozitioner). Его нужно перевернуть (Edit → Rotate 180).

1.2. Задать условия моделирования (Simulate → Simulation Setup): в окне Simulation Setup задать шаг Step Size – 0.0001, время моделирования Range End – 2.1. В области Integration Algoritm активировать метод Рунге-Кутта 4-го порядка. Запустить моделирование. Считаем, что полученное решение близко к точному. Замерить значения процесса при t = 1 с и 2 с (с точностью до четвертой значащей цифры после запятой). Для более точного измерения использовать “лупу”. Для этого надо выполнить следующие действия: нажать и удерживать клавишу Ctrl; с помощью указателя мыши выделить требуемую область осциллограммы; отпустить клавишу мыши; отпустить клавишу Ctrl. Для считывания координат использовать перекрестие (Правой клавишей мыши щелкнуть на Plot; в раскрывшемся окне левой кнопкой мыши щелкнуть на Read Coordinates). Если требуется возвратиться к исходному масштабу, проделать следующие действия: нажать и удерживать клавишу Ctrl; навести указатель мыши на окно; щелкнуть правой кнопкой мыши (или двойной щелчок левой); отпустить клавишу Ctrl.