Размещено на http://www.allbest.ru/

Министерство образования и науки Российской Федерации

Федеральное государственное автономное образовательное учреждение высшего профессионального образования

«Уральский федеральный университет имени первого Президента России Б.Н. Ельцина»

Институт радиоэлектроники и информационных технологий-РтФ

Кафедра Автоматики

ВЫПУСКНАЯ КВАЛИФИКАЦИОННАЯ РАБОТА БАКАЛАВРА

СТЕНД ЛАБОРАТОРНЫЙ: УПРАВЛЯЕМАЯ СВЯЗАННАЯ СИСТЕМА БАКОВ. МАТЕМАТИЧЕСКОЕ ОПИСАНИЕ. ИДЕНТИФИКАЦИЯ ПАРАМЕТРОВ МАТЕМАТИЧЕСКОЙ МОДЕЛИ

Руководитель, доц., к.т.н.

Е.Э. Страшинин

Студент гр. Р-400101

А.В. Хорошенин

Екатеринбург - 2014



Введение

В настоящей выпускной квалификационной работе, был подробно исследован и описан лабораторный стенд «Управляемая связанная система баков», спроектированный фирмой Feedback.

Исследовательская работа проводилась студентами гр. Р-400101 Хорошениным А.В. и Журавлевым А.П. совместно, в дальнейшем она была разделена на две части: первая часть - это данная работа, включает следующие разделы: описание стенда, математическое описание, идентификация параметров модели, проверка полученной модели; вторая часть: управление уровнями воды в одиночных и связанных баках.

В ходе анализа математической модели стенда, описанной в описании стенда фирмы Feedback[1], выяснилось, что в модели не учтены различные факторы (давление входящего напора жидкости, атмосферное давление, статическое давление жидкости и т.д.) и динамика системы описана с большими ограничениями (не учитывается турбулентное движение жидкости), что приводит к несоответствию математической модели и реального объекта.

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

Актуальность проблемы обусловлена тем, что проведенная работа может использоваться в качестве пособия для студентов по идентификации параметров объекта и составлению его математического описания, а также, как пособие для работы в MatLabSimulink.



1. Описание стенда

1.1 Описание установки

Лабораторный стенд представляет собой 4 цилиндрических сосуда (1), расположенных на вертикальной стойке (2). Каждый бак имеет 2 выходных отверстия (7) и (9). Отверстия (9) возможно перекрывать краниками, отверстия (7) - насадками, так же можно уменьшать площадь сечения отверстий (7) при помощи специальных насадков (см. рис. 2).К каждому сосуду подходят трубки (10) от насосов (4), расположенных в нижнем резервуаре для жидкости (3). Между верхними баками расположена трубка (11), при помощи которой реализуется система сообщающихся сосудов. Стенд включается при помощи кнопки питания (11). В каждом баке расположена сквозная металлическая трубка (6) высотой 27 см., выше данной трубки уровень воды в сосуде не поднимется. К каждому баку подключен датчик уровня (12), датчики установлены на обратной стороне стенда (рис. 2).

Рис. 1.1 - Схема стенда (передняя сторона)

С обратной стороны на стенде расположеныкалибруемые датчики уровня жидкости (12), модуль питания и усилитель мощности (14) PSUPA - PowerSupplyUnitandPowerAmplifier, так же данный модель преобразует управляющий сигнал и передаёт его на насосы, получает сигналы с датчиков уровня жидкости и передает их на ПК (15). Обмен информацией между стендом и ПК происходит через плату сопряжения (13) UCI (UniversalControlInterface). ПК соединяется с UCI по интерфейсу SCSI, UCI соединяется с PSUPA по интерфейсу LPT. Управляющие сигналы, которые являются напряжениями в диапазоне от 0 до 5В, передаются на усилитель мощности PSUPA, где они преобразуются в 24В управляющие сигналы. На управляющем ПК установлен пакет прикладных программ Matlab и среда Simulink.

Рис. 1.2 - Схема стенда (обратная сторона), соединение с ПК.

Рис. 1.3 - Насадок с маленьким отверстием. С площадью сечения отверстия s. Вид сбоку и вид сверху.

1.2 Программное обеспечение лабораторного стенда

Для обеспечения работы стенда и его взаимодействия с ПК необходим минимальный набор ПО:

- операционная система Windows XP;

- драйвер интерфейсной карты AdvantechPCI-1711;

- пакет прикладных программ Matlab 2009b с встроенной средой Simulink.

2. Математическое описание

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

2.1 Математическое описание одиночного сосуда

Таблица 2.1 - Перечень параметров объекта

Название параметра	Обозначение	Исходное значение	Единица измерения
Постоянная, связывающая управляющий сигнал и водяной поток	з	0.002076	м3/(В*с)
Площадь поперечного сечения выходной трубки без насадка, площадь поперечного сечения питающей трубки	а1	0.000050265	м2
Площадь поперечного сечения трубки с насадком	s	0.00001256	м2
Площадь поперечного сечения бака	А	0.013756	м2
Гравитационная постоянная	g	9.81	м/с2
Высота расположения питающей трубки относительно дна бака1	hТ1	0.3	м

Рассмотрим одиночный бак без насадка (рис.2.1), в дальнейшем составим уравнение для этого же бака с насадком.

Qвх- расход воды на входе сосуда, м3/с;

Qвых - расход воды на выходе из сосуда, м3/с;

h - уровень воды в сосуде, м;

hТ -высота питающей трубки относительно дна бака, м;

z1 - расстояние до центра массы воды в баке, м;

z2 - расстояние до центра массы воды в выходной трубке, м;

A - площадь поперечного сечения сосуда, м2;

a - площадь поперечного сечения выходной трубки, м2;

-скорость потока жидкости, поступающего в бак из насоса, м/с;

- вектор средней скорости истечения жидкости через сечение1-1(первая цифра в индексе означает номер бака, вторая - сечение), м/с;

- вектор средней скорости истечения жидкости через выпускное отверстие бака, м/с;

dl1 - перемещение жидкости в баке за время dt, м;

Рис.2.1. Одиночный бак

dl2 - перемещение жидкости в выходной трубке за время dt, м;

В любой момент времени при уровне h, объем воды в сосуде равен:

Если, то V - постоянный, h - постоянный, следовательно,

.

Если, то V - изменяется, h - изменяется, имеем скорость изменения объема воды в сосуде:

.

.

Из уравнения имеем:

.

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

Вывод уравнения Бернулли

Выделим в сосуде, содержащем идеальную (несжимаемую) жидкость, участки ограниченные сечениями 1-1 2-2 и 1'-1' 2'-2'. Полная энергия малого элемента жидкости равна:

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

Сверху на объем жидкости в баке действует сила , а снизу (минус, потому что сила направлена в противоположную сторону).

Здесь Р1и Р2 - статическое давление в точке пространства, где расположен центр массы рассматриваемого объема жидкости.

Статическое давление P (Па), действующее в покоящейся жидкости, складывается из внешнего равномерного давления на жидкость P0 (например, замеренное барометром атмосферное давление) и давления собственного веса жидкости (весового давления). Таким образом:

,

где z -расстояние от выбранного горизонтального уровня (нулевого уровня) до центра тяжести объёма жидкости (в нашем случаеz1 = ),

-равномерное внешнее давление на столб жидкости в сосуде.

В свою очередь формула для верхней границы следующая:

,

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

,

где - расстояние от нулевого уровня до трубки, - скорость потока жидкости, втекающего в бак, - плотность жидкости, - текущая высота столба жидкости.

Значение для нижней границыравно - давление среды, в которую вытекает жидкость. В нашем случае:

P02=Pа.

Пусть верхняя граница жидкости сдвинулась на dl1,а нижняя на dl2 за время dt.

Условие несжимаемости жидкости:

.

Объёмы, как видно, бесконечно малые, дифференциальные. Их самих можно рассматривать как дифференциалы объёма всего большого элемента.

При перемещении жидкости c расходом Qвых от сечения 1-1 к сечению 2'-2' её полная энергия во втором, конечном сечении всегда будет меньше, чем в начальном сечении. Это связано с тем, что в реальной жизни существуют силы, которые препятствуют движению жидкости. Учтем эти потери в кинетической энергии . Введем понятие коэффициента Кориолиса (б),учитывающий неравномерность распределения скоростей и равный отношению действительной кинетической энергии к кинетической энергии, определяемой по средней скорости. Также введем понятие коэффициента местного сопротивления отверстия(о), который учитывает внезапное изменение формы потока, скорости или направления ее движения. Тогда дифференциал кинетической энергии будет равен:

.

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



.

Учитывая формулу сократим на дифференциал объема:

.

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

.

Минус перед последним слагаемым получается из-за того что нулевой уровень находится выше, чем центр массы второго элемента

Распишем P1 и P2, используя формулы:

.

И приведем подобные:

.

На пути движения от начального сечения до конечного форма поперечных сечений потока может меняться самым причудливым образом. Однако, тот объем жидкости, который прошел за время t через любое сечение, должен остаться неизменным (в противном случае между сечениями изменялось бы количество проходящей жидкости, и течение было бы нестационарным). Стационарное течение жидкости - это такое течение, при котором скорость жидкости в каждой данной точке остается постоянной как по величине, так и по направлению. Для такого режима форма и расположение линий течения со временем не изменяются. Таким образом:

,

,

Подставим полученное выражение скорости в уравнение и выразим скорость истечения жидкости из бака (:

,

,

=,

,

,



.

.

В уравнении величинаназывается коэффициентом по скорости. Таким образом, скорость истечения жидкости из бака, равна:

.

Уравнение приняло вид закона Торричелли. На основании уравнения запишем расход на выходе:

,

,

где м1 = 1.33 - коэффициент расхода отверстия первого бака равен отношению коэффициента расхода насадка и коэффициента расхода отверстия в дне бака.

Так как насадки и отверстия у баков идентичные, то (.

считается по следующей формуле:

,

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

Окончательно уравнение примет вид:

Уравнение отражает скорость изменения уровня воды h в одиночном баке без насадка.

Для составления уравнения, характеризующего скорость изменения воды в баке с насадком, необходимо в уравнении изменить диаметр отверстия трубки на диаметр отверстия насадка. Тогда уравнение примет следующий вид:

2.2 Математическое описание двух связанных баков

Для составления математического описания двух связанных баков - сообщающихся сосудов была рассмотрена система, представленная на рисунке 2.2. Кран, перекрывающий перешеек между баками 1 и 3, открыт, значит, имеется влияние уровня жидкости в одном баке на уровень жидкости в другом и необходимо учитывать переток жидкости через перешеек. Так же уравнения были составлены и для случая, когда на выходных трубках находился насадок.



Рис. 2.2 Система связанных баков

Перечень параметров системы приведен в таблице 2.2

Таблица 2.2 - Перечень параметров объекта

Название параметра	Обозначение	Исходное значение	Единица измерения
Постоянная, связывающая управляющий сигнал и водяной поток	з	0.002076	м3/(В*с)
Площадь поперечного сечения выходной трубки из бака 1	а1	0.000050265	м2
Площадь поперечного сечения выходной трубки из бака 3	а2	0.000050265	м2
Площадь поперечного сечения выходной трубки из бака 1	А1	0.000050265	м2
Площадь поперечного сечения выходной трубки из бака 3	А3	0.000050265	м2
Площадь перешейка	Sп	0.000050265	м2
Площадь поперечного сечения насадка	s	0.00001256	м2
Площадь поперечного сечения баков	А	0.013756	м2
Гравитационная постоянная	g	9.81	м/с2
Высота расположения питающей трубки относительно дна бака1	hТ1	0.3	м
Высота расположения питающей трубки относительно дна бака 3	hТ3	0.3	м



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

,

.

Для определения количества жидкости, перетекаемой через перешеек необходимо ввести коэффициент перетока б, а так же площадь перешейка Sп.

Тогда система уравнений, описывающая сообщающиеся сосуды примет вид 2.25:

,

.

Знак «» отражает направление перетока, т.е. жидкость перетекает из бака с большим уровнем в бак с меньшим. Баку, уровень жидкости в котором выше, соответствует знак «-», то есть жидкость из него вытекает. Баку, уровень жидкости в котором ниже, соответствует знак «+», то есть жидкость в него втекает через перешеек.

Для составления уравнений, характерных системы с насадками, необходимо изменить площадь выходных отверстий на площадь отверстий насадков. Тогда система 2.25 примет следующий вид:



,

.



3. Определение статической характеристики двигателей

Рассмотрим систему уравнений (2.24), полученную в предыдущем разделе. Переменные, входящие в нее, описываются через поток жидкости, входящий в первый и третий баки - и . Найдем зависимость скорости потока жидкости от напора (статическую характеристику двигателя).

Статическая характеристика элемента - это зависимость установившихся значений выходной величины от значения на входе системы. Входной величиной для двигателя является напор Qвх, а выходной - скорость потока жидкости.

Скорость потока жидкости определяется по следующей формуле:

,

Где t - время наполнения бака до уровня h, - объем бака.

Для определения статической характеристики двигателей был проведен ряд экспериментов в модели реального времени (рис. 3.1.) с первым и вторым насосами. В проведенных экспериментах определялось время наполнения баков до уровня 20 см на различных напряжениях: U = 2.5, 2.75, 3, 3.25, 3.5 В. В результате проведенных экспериментов были получены зависимости скорости наполнения бака от входного напора (таблицы 3. 1 и 3.2).

Рис. 3.1 Эксперимент по измерению скорости наполнения первого бака



Таблица 3.1 - Зависимость скорости наполнения первого бака от входного напора.

U, В	t1, с	V1, м/с
2.5	39.7	1.378688
2.75	33.5	1.633848
3	29.5	1.855387
3.25	26.5	2.0654306
3.5	22.5	2.432618

t1 - время наполнения первого бака до уровня h = 0.2 м в зависимости от напряжения, подаваемого на насос U, V1 - скорость входного потока.

Таблица 3.2 - Зависимость скорости наполнения второго бака от входного напора.

U, В	t1, с	V1, м/с
2.5	43.3	1.264063
2.75	36.8	1.487335
3	32	1.710435
3.25	28.5	1.920488
3.5	26.2	2.089081

t2 - время наполнения первого бака до уровня h = 0.2 м в зависимости от напряжения, подаваемого на насос U, V2 - скорость входного потока.

По полученным результатам была проведена обработка результатов характеристики с помощью инструмента в Matlab Curve Fitting Tool(cftool).

В результате были определены законы изменения скорости вытекания жидкости из питающей с применением линейных полиномов для первого насоса и для второго (рис. 3.2 и 3.3).

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

() = 1.016 - 1.174



где u(t) - напряжение подаваемое на первый насос, () - скорость потока жидкости, питающее первый бак.

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

() = 0.8333 - 0.8055,

где u(t) - напряжение подаваемое на второй насос, () - скорость потока жидкости, питающей второй бак.

Рис. 3.2 - Статическая характеристика насоса (бак 1)

Рис. 3.3 - Статическая характеристика насоса (бак 3)

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



4. Идентификация параметров нелинейной модели

Для проведения идентификации параметров нелинейной модели была использована методика, разработанная магистром кафедры Автоматики Чебыкиным Д.В.

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

Этот раздел посвящён описанию некоторых методов решения задач идентификации с помощью графического интерфейса System Identification Tool GUI пакета MATLAB.

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

Таким образом, встает задача определения параметров модели. Один способ их нахождения - это непосредственное измерение их величин на физическом объекте или использование стандартных параметров.

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

4.1 Организация массива данных

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

Рис.4.1 - Постановка эксперимента для идентификации

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

В этом эксперименте на вход объекта подается постоянное значение напряжения. Это проводится с целью упрощения процесса идентификации, таким образом, проводится идентификация параметров в статике и на основе этих данных строятся зависимости для всех значений управляющего сигнала. Реакция объекта на этот сигнал должна быть записана и сохранена в рабочей памяти MATLAB. Это осуществляется с помощью модуля ToWorkspace. В параметрах этого модуля должны быть указаны его имя, интервал дискретизации по времени и формат записи (удобно - Array (рис. 4.2)). Следует продуманно выбрать интервал дискретизации Т0. Его величина должна быть достаточно мала, чтобы записанный процесс отражал всю суть протекающего процесса. Одновременно с этим он должен быть достаточно велик, чтобы не записывать помехи сигнала. Это может привести впоследствии к грубости синтезированной системы управления. Интервал дискретизации задаётся в модуле ToWorkspace свойством SampleTime.

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

По окончании эксперимента в массиве Simout будет сформировано три столбца данных - управление u, уровень воды в баке 1 - h1 и уровень воды в баке 2 - h2.

Массив Simout следует сохранить в виде mat-файла: установив курсор на значке Simout в окне Workspace, после нажатия правой клавишей мыши выбрать `SaveAs' и написать имя mat-файла, например, `Experiment'.

Рис.4.2 - Параметры модуля ToWorkspace

4.2 Идентификация

Модель исследуемого реального объекта описывается набором (вектором) параметров. Для того чтобы данная модель в достаточной степени была адекватна реальному объекту, значения этих параметров должны соответствовать действительным параметрам физического объекта. Таким образом, встаёт задача определения параметров модели. Один способ их нахождения - это непосредственное измерение их величин на физическом объекте или использование стандартных (справочных) параметров.

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

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

Таким образом, с математической точки зрения, идентификация параметров - это задача оптимизации. Пусть u(t)- пробное воздействие (или управление), которое подаётся на реальный управляемый физический объект, а y(y) - реакция объекта, измеренная на его выходе. Заметим, что значения процесса y(t) являются экспериментальными данными, которые несут в себе информацию о конкретном поведении управляемого объекта.

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

где [0,t] - интервал времени, на котором осуществлялось измерение выхода управляемого объекта.

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

В пакете Simulink имеется встроенное средство (процедура) Simulink Design Optimization для идентификации параметров модели. Эта процедура проста и реализует процесс минимизации и является достаточно прозрачной и понятной для пользователя.

Процедура Simulink Design Optimization поддерживает различные источники данных для идентификации моделей. В том числе, зависимости u(t) и могут быть получены прямо с реального оборудования. Для этого необходимо поставить эксперимент, который будет возбуждать динамику системы, зависящую от идентифицируемых параметров. Этот процесс, представляющий собой выход установки, необходимо записать в рабочую область Matlab с помощью блока ToWorkspace, который находится в библиотеке Simulink (SimulinkLibraryBrowser/Sinks/ToWorkspace).

В свойствах данного блока (Рис. 4.2.) в поле SaveFormat: выбрать из выпадающего списка пункт Array, в окне SampleTime записать шаг дискретизации (на рисунке - 0,1 с). На вход этого блока через модуль объединения нескольких сигналов в вектор Mux (SimulinkLibraryBrowser/SignalRouting/Mux) подать сигнал управления u и выходные сигналы объекта (например, h1 и h2) на модуль ToWorkspace, как показано на рис. 4.3.

Рис.4.3 - Подключение сигналов к блоку ToWorkspace

После окончания эксперимента нужно сохранить полученный массив данных. Для этого в окне WorkspaceMatlab навести курсор на иконку Simout и сохранить как mat-файл, выбрав имя, например, “Experiment”. В оболочке Simulink создать рабочий проект идентификации. В рабочее окно Matlab поместить файл “Experiment.mat”. Там же должна располагаться модель, параметры которой будут идентифицироваться. Все параметры этой модели, должны быть определены в m-файле “Modelfile”. В том числе в этом файле должны быть указаны идентифицируемые параметры, для которых должны быть указаны их предварительные ориентировочные значения. В файл должны быть включены операторы загрузки файла “Experiment.mat”:

load-matExperiment%Протокол эксперимента на стенде

и чтения из него записанных в эксперименте сигналов:

u = simout(:,1);

h1 = simout (:,2);

h2 = simout (:,3);

Этот файл должен быть прикреплён к mdl-файлу модели. Для этого его имя следует вписать (без расширения) в окно Modelpre-loadfunction (File/ModelProperties/Callbacks/PreloadFcn).

Перед началом процесса идентификации параметров, необходимо:

· Отключить управляющие сигналы от входа модели. Вместо них подключить к входам модели (если их несколько) специальные компоненты, представляющие входные порты, которые находятся в библиотеке Simulink (SimulinkLibraryBrowser/CommonlyUsedBlocks/In или Simulink/Sinks/In).

· К выходу модели подключить соответствующие компоненты, представляющие выходные порты, которые находятся также в библиотеке (SimulinkLibraryBrowser/CommonlyUsedBlocks/Out или Simulink/Sinks/Out).

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

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

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

4.3 Порядок проведения идентификации

Чтобы начать процесс идентификации параметров, необходимо открыть mdl-файл модели и из меню Simulink: Tools/Parameter Estimation запустить опцию Controland Estimation Tools Manager. В результате появится окно процедуры идентификации параметров, показанное на рис.4.4.

Далее, щелчком в дереве (рис. 3.4.) слева по кнопке TransientData открываем окно Transientdatasets и нажимаем кнопку New. В результате в окне Propertiesпоявляется имя рабочего сеанса «Newdata», которое по желанию можно заменить на любое другое.

После этого, нажимая слева в дереве на кнопку Newdata, вводим в правое окно данные измерений: в соответствии с предыдущим примером, на вкладке Inputdata - в соответствующие столбцы имя управляющего сигнала u и период дискретизации (0.1); на вкладке Ouputdata -h1, h2ипериоды дискретизации (0.1, 0.1).

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

Рис.4.4 - Окно Control and Estimation Tools Manager

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

Нажав на кнопку Pre-process… (Рис. 3.5.) можно настроить введенные данные. Перейдя на вкладку ExclusionRules (Правила исключения) можно выбрать один или несколько из следующих пунктов для корректировки значений процесса:

Bounds - можно выбрать верхние и нижние границы для времени и значений данных.

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

Matlabexpression - логическое выражение Matlab, которое выбирает конкретные значения данных.

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

Рис.4.5 - Окно Data Preprocessing Tool for Dataset: Exclusion Rules

НавкладкеDetrend/Filtering (Рис.3.6.)имеются следующие возможности:

Поставив галочку на Detrending, вы можете выбрать постоянное или прямое удаление тренда. Постоянное (Constant) удаление тренда удаляет среднее значение данных для создания нулевого среднего данных. Прямая линия (StraightLine) удаления тренда находит линейные тренды (методом наименьших квадратов) и затем удаляет их.

Выбрав пункт Filtering, у вас появляется возможность отфильтровать ваши данные следующим образом:

FirstOrder - фильтр типа ,где ф - постоянная времени, указываемая в соответствующем поле;

TransferFunction - фильтртипа, где вы задаёте коэффициенты в виде векторов в соответствующих полях коэффициентов A и B;

Ideal - идеализированный вариант фильтра, где устанавливается либо полоса пропускания (Passband), либо полоса затухания (Stopband). Указывается фильтр из двух элементов вектора в поле Range (Hz) (Диапазон). Эти фильтры являются идеальными в том смысле, что не существует конечного спада частотной характеристики или колебаний, т.е. концы диапазона горизонтальны в частотнойобласти.

Рис.4.6 - Окно Data Preprocessing Tool for Dataset: Detrend/Filtering

Нажав на кнопку ExcludeGraphically можно откорректировать процесс,используя графический интерфейс.

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

После корректировки значений необходимо выбрать, куда сохранять откорректированные данные вверху окна (Рис.3.6.). ВполеWriteResultsto: выбираем один из пунктов existingdataset, что заменит набор данных на измененный, или newdataset, что создаст новый набор данных с введенным именем. По окончании корректировки каждого из сигналов следует восстанавливать имя изменённого набора данных (например, Dataset1) и нажимать кнопку Add. Аналогично вводятся и корректируются данные измерений, которые соответствуют выходу модели. После окончания корректировки последнего сигнала закрыть это окно.

ТеперьвлевомокненавкладкеControlandEstimationToolsManagerвдереве TransientData появится две вкладки: NewData - с исходными данными и Dataset1 - откорректированными данными (рис.3.7.).

Рис.4.7 - После коррекции сигналов

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

- в дереве операций (Рис.3.7.) необходимо выделить пункт Variables;

- далее справа нажать кнопку Add…;

- в открывшемся окне SelectParameters (рис.3.8.) выделять переменные, подлежащие идентификации из тех, которые были объявлены в файле”Modelfile.m”;

- для каждой выделенной переменной нажимать кнопку Applyили выделить все сразу, зажав клавишу “Ctrl”;

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

- после введения имени последней переменной нажать кнопку ОК.

В данном окне нажав на определенный параметр можно просмотретьего значение, установить Initial guess (начальное значение), интервалы изменения(Minimum, Maximum), Typical value (типичное значение).

Рис.4.8 - Задание идентифицируемых параметров модели

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

- в дереве (Рис.3.7.) выделить пункт Estimation;

- нажать справа кнопку New, создав тем самым новую задачу идентификации параметров с именем по умолчанию NewEstimation;

- выделив пункт NewEstimation в дереве, под пунктом Estimation справа получим несколько вкладок, в которых настраивается и осуществляется процесс идентификации параметров.

На вкладке DataSets (Рис.3.9.) необходимо установить флажок в списке наборов данных, т.е. выбрать из этого списка созданный ранее набор данных NewData, например, Dataset1. При этом для идентификации параметров будет использоваться именно этот набор данных.



Рис.4.9 - Выбор данных для идентификации параметров

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

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

Рис.4.10 - 0Выбор идентифицируемых параметров

После выполнения всех операций предварительной настройки, можно переходить на последнюю вкладку Estimation.

Нажав на кнопку EstimationOptions, можно настроить разные параметры процесса идентификации:

- время моделирования процесса, по которому происходит идентификация;

- величина временного шага;

- методы и алгоритмы оптимизации (подробнее по методам и алгоритмам оптимизации см. ProductHelp MATLAB).

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

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

Процедура идентификации завершаетсяпри:

- достижении функционалом минимума с заданной точностью;

- достижении минимального заданного шага изменения идентифицируемых параметров;

- достижении заданного числа итераций;

- нажатии на кнопку Stop(принудительно).

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

Значения компонент вектора идентифицированных параметров модели автоматически заносятся в рабочую область. Кроме того, эти значения можно посмотреть на вкладке Parameters.



Рис.4.11 - Результаты процесса идентификации

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

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

a. в дереве (рис.3.9.) нужно выбрать пункт Estimation/Views;

b. нажать на кнопку New; при этом создается некоторое количество графиков процессов;

c. выделив созданный пункт NewViews в дереве, справа получить таблицу-перечень графиков (Рис.3.12.), доступных для визуализации;

d. далее из выпадающего меню можно выбрать конкретные графики, которые программа должна визуализировать.

Возможен вывод следующих графиков, иллюстрирующих протекание процесса идентификации:

- Measuredandsimulated - два наложенных графика - измеренный реальный процесс и модельный процесс на выходе модели при найденном векторе её параметров;



Рис.4.12 - Визуализация решения

- Costfunction - показывает погрешность невязки этих двух процессов, т.е. минимальное полученное значение критерия ;

- Parametersensitivity - представляет чувствительность модели по идентифицированным параметрам;

- Parametertrajectory - показывает траекторию изменения параметров в процессе идентификации;

- Residuals - разность между измеренным реальным процессом и модельным процессом на выходе модели при найденном векторе оценок параметров объекта.

4.4 Подготовка к проведению идентификации

Следует отметить, что идентификация параметров объекта проводилась при наличии насадков на вхдных трубках и без насадков. Обусловлено это тем, что параметры таких систем различные.

Также для дальнейшей идентификации параметров модели, необходимо её упростить.

Из уравнений 2.22 и 2.23 видно, что идентифицировать придется массив параметров: . Для уменьшения количества перемеренных обозначим произведение коэффициентов одной переменной: и , где - общий коэффициент по первому баку, - общий коэффициент по третьему баку.

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

Рис. 4.13 - Общая модель системы

Уровни воды в баках до и после упрощения полностью совпадают.

Составим упрощенную модель:

,

.

4.5 Результаты идентификации

В нелинейной модели необходимо идентифицировать четыре параметра, а именно: отдельно для каждого бака с насадком на трубке и без. Для этого воспользуемся инструментом SystemIdentificationTool GUI пакета Matlab, подробное описание работы с которым, приведено выше.

Первый бак, трубка без насадка

Для процедуры идентификации нами был проведен ряд экспериментов для получения данных с датчиков уровня жидкости при напряжениях 2.5; 2.7; 3; 3.3; 3.5; 3.8; 4; 4,6; 5 В. При напряжениях меньше 2.5 В, уровень воды в баке минимальный.

Продолжительность измерений соответствовала значению до 500с (5000 тактов), так как за это время при каждом из напряжений ниже 3.8 В наблюдался установившийся режим. Начиная с 3.8 В и до 5 В уровень жидкости в баке превышал критическую отметку в 27 см, и вода переливалась в сквозную трубку. Поэтому при высоких напряжениях эксперименты заканчивались в момент преодоления водой уровня 27 см.

Графики, полученные в результате данных измерений для первого бака, приведены на рисунках 4.14-4.31.

В результате экспериментов в пакете Matlab записывались следующие параметры (способ записи параметров описан в разделе 4.1): время, уровень воды в первом баке, напряжение, подаваемое на насос.

На основе полученных данных инструмент SystemIdentificationTool GUI может подобрать такие параметры которые будут максимально точно отображать модель.

Полученные данные в виде таблицы будут приведены в разделе 4.7.

Рис. 4.14 - Результат идентификации параметров при напряжении, подаваемом на насос u=2.5 В



Рис. 4.15 - Параметрическая траектория

Рис. 4.16 - Результат идентификации параметров при напряжении, подаваемом на насос u=2.7 В

Рис. 4.17 - Параметрическая траектория

Рис. 4.18 - Результат идентификации параметров при напряжении, подаваемом на насос u=3 В

Рис. 4.19 - Параметрическая траектория

Рис. 4.20 - Результат идентификации параметров при напряжении, подаваемом на насос u=3.3 В

Рис. 4.21 - Параметрическая траектория

Рис. 4.22 Результат идентификации параметров при напряжении, подаваемом на насос u=3.5 В

Рис. 4.23 - Параметрическая траектория

Рис. 4.24 - Результат идентификации параметров при напряжении, подаваемом на насос u=3.8 В

Рис. 4.25 - Параметрическая траектория

Рис. 4.26 - Результат идентификации параметров при напряжении, подаваемом на насос u=4 В



Рис. 4.27 - Параметрическая траектория

Рис. 4.28 - Результат идентификации параметров при напряжении, подаваемом на насос u=4.6 В

Рис. 4.29 - Параметрическая траектория

Рис. 4.30- Результат идентификации параметров при напряжении, подаваемом на насос u=5 В

Рис. 4.31 - Параметрическая траектория

Первый бак, трубка с насадком

Для процедуры идентификации нами был проведен ряд экспериментов для получения данных с датчиков уровня жидкости при напряжениях 2.1; 2.5; 3; 3.5; 4; 5 В. Продолжительность измерений соответствовала значению времени, при котором уровень воду в баке равнялся критическому (27 см).

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

В результате экспериментов в пакете Matlab записывались следующие параметры (способ записи параметров описан в разделе 4.1): время, уровень воды в первом баке, напряжение, подаваемое на насос.

На основе полученных данных инструмент SystemIdentificationTool GUI может подобрать такие параметры

Полученные данные в виде таблицы будут приведены в разделе 4.7.

Рис. 4.32 Результат идентификации параметров при напряжении, подаваемом на насос u=2.1 В



Рис. 4.33 - Параметрическая траектория

Рис. 4.34 - Результат идентификации параметров при напряжении, подаваемом на насос u=2.5 В

Рис. 4.35 - Параметрическая траектория

Рис. 4.36 Результат идентификации параметров при напряжении, подаваемом на насос u=3 В

Рис. 4.37 - Параметрическая траектория

Рис. 4.38 - Результат идентификации параметров при напряжении, подаваемом на насос u=3.5 В

Рис. 4.39 - Параметрическая траектория

Рис. 4.40 - Результат идентификации параметров при напряжении, подаваемом на насос u=4 В



Рис. 4.41 - Параметрическая траектория

Рис. 4.42 - Результат идентификации параметров при напряжении, подаваемом на насос u=5 В

Рис. 4.43 - Параметрическая траектория

Третий бак, трубка без насадка

Для процедуры идентификации был проведен ряд экспериментов для получения данных с датчиков уровня жидкости при напряжениях 2.5; 2.7; 3; 3.3; 3.5; 3.8; 4; 4,6; 5 В. При напряжениях меньше 2.5 В, уровень воды в баке минимальный.

Продолжительность измерений соответствовала значению до 500с (5000 тактов), так как за это время при каждом из напряжений ниже 3.8 В наблюдался установившийся режим. Начиная с 3.8 В и до 5 В уровень жидкости в баке превышал критическую отметку в 27 см, и вода переливалась в сквозную трубку. Поэтому при высоких напряжениях эксперименты заканчивались в момент преодоления водой уровня 27 см.

В результате экспериментов в пакете Matlab записывались следующие параметры (способ записи параметров описан в разделе 4.1): время, уровень воды в первом баке, напряжение, подаваемое на насос.

На основе полученных данных инструмент SystemIdentificationTool GUI может подобрать такие параметры которые будут максимально точно отображать модель.

Полученные данные в виде таблицы будут приведены в разделе 4.7.

Рис. 4.44 - Результат идентификации параметров при напряжении, подаваемом на насос u=2.5 В

Рис. 4.45 - Параметрическая траектория

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

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

Рис. 4.46 - Результат идентификации параметров при напряжении, подаваемом на насос u=2.7 В

Рис. 4.47 - Параметрическая траектория

Рис. 4.48 - Результат идентификации параметров при напряжении, подаваемом на насос u=3 В



Рис. 4.49 Параметрическая траектория

Рис. 4.50 - Результат идентификации параметров при напряжении, подаваемом на насос u=3.3 В

Рис. 4.51 - Параметрическая траектория

Рис. 4.52 - Результат идентификации параметров при напряжении, подаваемом на насос u=3.5 В



Рис. 4.53 - Параметрическая траектория

Рис. 4.54 - Результат идентификации параметров при напряжении, подаваемом на насос u=3.8 В

Рис. 4.55 - Параметрическая траектория

Рис. 4.56 Результат идентификации параметров при напряжении, подаваемом на насос u=4 В

Рис. 4.57 - Параметрическая траектория

Рис. 4.58 - Результат идентификации параметров при напряжении, подаваемом на насос u=4.6 В

Рис. 4.59 - Параметрическая траектория

Рис. 4.60 Результат идентификации параметров при напряжении, подаваемом на насос u=5 В

Рис. 4.61 - Параметрическая траектория

Третий бак, трубка с насадком

Для процедуры идентификации нами был проведен ряд экспериментов для получения данных с датчиков уровня жидкости при напряжениях 2.1; 2.5; 3; 3.5; 4; 5 В. Продолжительность измерений соответствовала значению времени, при котором уровень воду в баке равнялся критическому (27 см).

В результате экспериментов в пакете Matlab записывались следующие параметры (способ записи параметров описан в разделе 4.1): время, уровень воды в первом баке, напряжение, подаваемое на насос.

На основе полученных данных инструмент SystemIdentificationTool GUI может подобрать такие параметры

Полученные данные в виде таблицы будут приведены в разделе 4.7.

Рис. 4.62 - Результат идентификации параметров при напряжении, подаваемом на насос u=2.1 В



Рис. 4.63 - Параметрическая траектория

Рис. 4.64 - Результат идентификации параметров при напряжении, подаваемом на насос u=2.5 В

Рис. 4.65 - Параметрическая траектория

Рис. 4.66 - Результат идентификации параметров при напряжении, подаваемом на насос u=3 В

Рис. 4.67 - Параметрическая траектория

Рис. 4.68 - Результат идентификации параметров при напряжении, подаваемом на насос u=3.5 В

Рис. 4.69 - Параметрическая траектория

Рис. 4.70 - Результат идентификации параметров при напряжении, подаваемом на насос u=4 В



Рис. 4.71 - Параметрическая траектория

Рис. 4.72- Результат идентификации параметров при напряжении, подаваемом на насос u=5 В

Рис. 4.73 - Параметрическая траектория

Коэффициент перетока

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

4.6 Обработка данных идентификации

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

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

Таблица 4.1 - Результаты идентификации для 1 бака без насадков

Параметр/ напряжение			Невязка
2.5	2	0.01287500	50
2.7	1.16540	0.12795000	430
3.0	1.02890	0.15710000	240
3.3	1.00210	0.15489000	599
3.5	0.98559	0.15461000	833
3.8	1.05160	0.08928500	254
4.0	1.13630	0.01851800	163
4,6	1.23740	0.00031616	45
5.0	1.27410	0.00078261	24

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

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

Результат аппроксимации параметра представлен на рисунке 4.74. Уравнение, соответствующее зависимости параметра от входного напора:

Рис. 4.74 - Зависимость от Qвх1

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

Рис. 4.75. Зависимость от Qвх1 (полином второй степени)



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

Для параметра также получен полином третьей степени.

График, полученный в результате аппроксимации параметра , представлен на рисунке 4.76.

Для бака, на выходной трубке которого присутствует насадок, получены значения, приведенные в таблице 4.2.

f(x) = p1*x^3 + p2*x^2 + p3*x + p4

Coefficients (with 95% confidence bounds):

p1 = 1.857e+012 (-5.486e+011, 4.263e+012)

p2 = -6.758e+008 (-1.542e+009, 1.907e+008)

p3 = 8.054e+004 (-1.968e+004, 1.808e+005)

p4 = -2.149 (-5.868, 1.571)



Linear model Poly3:

f(x) = p1*x^3 + p2*x^2 + p3*x + p4

Coefficients (with 95% confidence bounds):

p1 = -1.313e+012 (-3.07e+012, 4.444e+011)

p2 = 5.105e+008 (-1.224e+008, 1.143e+009)

...