Автоматизированное проектирование нелинейных систем управления

[Введите текст]

Задание

1. Построение модели диагностического стенда Вычислить по паспортным данным двигателя необходимые параметры электромеханической системы «двигатель?генератор» и построить его модель в среде Simulink. Паспортные данные электродвигателя приведены в таблице.

Таблица 1 — Технические данные тягового двигателя постоянного тока

2. Получение модели устройства управления По передаточной функции регулятора, используя аппарат переменных состояния, получить модель «вход? состояние? выход». Также привести схему имитационного моделирования ПИД-регулятора в Simulink.

3. Определение параметров регулятора Построить в Simulink схему имитационного моделирования замкнутой системы и осуществить настройку (определение численных значений параметров) регулятора по минимуму динамической ошибки в системе. Промоделировать работу системы и построить временные зависимости для управляющих напряжений u1 и u2 и выходной переменной y (t)? скорости вращения (t).

4. Построение модели замкнутой системы Записать математические модели ПИД-регулятора и объекта относительно входного сигнала g (t) и вектора состояния x (t). Получить модель «вход? состояние? выход» для расширенного вектора состояния z (t) = |x (t) v (t)|? замкнутой системы управления.

5. Имитационное моделирование системы управления Составить в Matlab программы вычисления правых частей обыкновенных дифференциальных уравнений и численного интегрирования методом Рунге-Кутта четвертого порядка. Выбрать шаг интегрирования h из условия обеспечения численной устойчивости решения. Составить программу для сравнения результатов численного решения дифференциальных уравнений в соответствии с разработанной программой и стандартной процедурой ode45. При выбранном шаге h провести численное решение и сравнить результаты моделирования.

Для выходной переменной y (t) сравнить с результатами имитационного моделирования замкнутой системы в Simulink.

6. Проведение статистического анализа Построить в Simulink схему имитационного моделирования замкнутой системы со случайным возмущением x (t) = Mвозм (t), используя блок генерации независимых нормально распределенных случайных сигналов. Провести имитационное моделирование и статический анализ случайных величин. Получить статистические характеристики для сигналов x (t) и y (t). Сравните оценки математического ожидания и дисперсии сигнала x (t) с заданными параметрами блока генерации. Для данного числа интервалов группирования получите и построить эмпирические законы распределения (гистограммы). Провести имитационное моделирование и статистический анализ случайных процессов во временной области. Получить и построить оценки автокорреляционных функций bRx (ф) и bRy (ф). Сравнить полученную bRx (ф) с идеальной характеристикой белого (независимого шума). Провести имитационное моделирование и статистический анализ спектральных характеристик случайных процессов. Сравнить выборочные оценки bSx (щ) с моделью белого шума.

7. Оформление и содержание пояснительной записки Пояснительная записка оформляется в полном соответствии со стандартами предприятия (университета).

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

Реферат

Автоматизированная система управления, двигатель, генератор, регулятор, устойчивость, скорость вращения, ток, напряжение, объект управления, устройство управления.

Объектом разработки является регулятор для диагностического стенда испытания тяговых двигателей.

Пояснительная записка выполнена в программе Microsoft Word 2007, расчеты выполнены в среде Matlab, моделирование произведено в пакете Simulink.

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

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

Математическое и программное обеспечение конкретных САПР в большей степени зависит от объекта проектирования.

При автоматизации проектирования САУ, где основным математическим описанием систем во временной области являются переменные состояния, записанные в векторной или матричной форме, широко используется Matlab .

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

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

Методы проектирования всегда являются многовариантными, особенно, при проведении статистического анализа.

1. Построение модели объекта проектирования

1.1 Описание объекта управления и его функциональных моделей

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

Две одинаковые электрические машины соединяются между собой механически и электрически и подключаются к внешнему источнику энергии. Одна из машин работает в режиме генератора и отдает всю вырабатываемую электрическую энергию другой машине, которая работает в режиме двигателя и расходует всю свою механическую энергию на вращение первой машины. Расход энергии при испытаниях по методу взаимной нагрузки определяется суммарными потерями в обеих машинах. Если учесть, что КПД тяговых электродвигателей превышает 90%, то оказывается, что для испытаний требуется источник мощности, составляющий всего 10? 20% мощности каждой испытуемой машины, в этом и заключается экономичность метода.

Рисунок 1 — Структурная (а) и принципиальная (б) схемы испытательной станции Для введения энергии в систему применяется способ параллельного включения источника, когда якорные обмотки машин включаются параллельно и когда подключаемый к ним линейный преобразователь (ЛП) обеспечивает необходимый режим напряжения. Компенсацию электрических потерь выполняет вольтодобавочный преобразователь (ВДП) путем регулирования тока в контуре «двигатель? генератор». Схемы испытательной станции приведены на рисунке 1.

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

1.2 Динамическая модель объекта

Электрическая схема замещения системы «двигатель? генератор» представлена на рисунке 2.

Рисунок 2 — Электрическая схема замещения В цепи двигателя включены последовательно обе обмотки двигателя (Д): якорная с сопротивлением Rя. д и индуктивностью Lя. д и обмотка возбуждения с сопротивлением Rв. д и индуктивностью Lв. д, а также обмотка возбуждения генератора (Г) с сопротивлением Rв. г и индуктивностью Lв.г. В цепи генератора включена только якорная обмотка генератора с сопротивлением Rя. г и индуктивностью Lя.г. ЭДС eд и eг действуют в якорных цепях двигателя и генератора, в которых также протекают токи iд и iг, соответственно. На основании второго закона Кирхгофа, записанного для контура К1, выполняется уравнение электрического баланса:

.

Уравнение электрического баланса для контура К2 записывается аналогично:

.

В уравнения (1) и (2) входят следующие величины: напряжения линейного u1 и вольтодобавочного u2 преобразователей, В; ЭДС двигателя eд и генератора eг, В; ток цепи двигателя iд, А; сопротивления Rв. д, Rя. д, Rв. г, Ом; индуктивности Lв. д, Lя. д, Lв. г, Гн.

Уравнение механического баланса получается на основании второго закона Ньютона и имеет вид:

где Щ? угловая скорость вращения валов электромеханической системы, рад/с; Mд? момент двигателя, Н· м; Mг? момент генератора, Н· м; Mв? суммарный механический момент внешних сил, действующих на вал двигателя и генератора, Н· м; J? момент инерции системы, кг· м2.

Введем обозначения:

; (4)

;(5)

;

и перепишем уравнения (1) — (3) для изображений сигналов:

;(8)

;(9)

.(10)

Значения вращающих моментов двигателя и генератора зависят от токов в их якорных обмотках и рассчитываются по формулам:

;(11)

а связь ЭДС двигателя и генератора с угловой скоростью вращения вала описывается соотношениями:

;

.(14)

Коэффициенты kM и kE зависят от конструктивных параметров электрических машин и тока iд, который протекает в их обмотках возбуждения. Качественный вид зависимостей этих коэффициентов от тока iд показан на рисунке 3 и может быть описан нелинейной функцией fв (iд).

Рисунок 3 — Качественный вид кривой намагничивания В курсовой работе такие нелинейные функциональные зависимости будем задавать функцией гиперболического тангенса:

(15)

где Imax? ток насыщения; б? параметр нелинейности, б = 2.

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

Ток насыщения примем равным максимальному току:

где Iн? номинальный ток двигателя, определяемый по формуле

где КПД? коэффициент полезного действия двигателя.

;

.

Запишем выражения для функций kM и kE:

;

.

где cM, cE? постоянные коэффициенты.

Тогда формулы (11) — (14) примут вид:

;(20)

;(21)

.

Значения коэффициентов cM и cE определим по паспортным данным двигателя. Рассмотрим схему с последовательным возбуждением ДПТ, которая приведена на рисунке 4.

Рисунок 4 — Схема включения двигателя с последовательным возбуждением Уравнение электрического баланса в установившемся (номинальном статическом) режиме будет иметь вид:

.(23)

Значение ЭДС может быть найдено при номинальной скорости вращения Щн по формуле:

(24)

тогда из выражений (23) и (24) следует:

;

.

В номинальном режиме связь вращающего момента с током описывается соотношением:

(26)

с другой стороны, номинальный момент определяется по выражению:

.

Из формул (26) и (27) следует, что

;

.

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

где p? число пар полюсов электродвигателя,

.

Для простоты будем считать, что все индуктивности равны: Lя. д = Lв. д == Lя. г = Lв.г.

Момент инерции выберем из диапазона:

;

;

;

Выберем значение момента инерции J = 4710 кгм2.

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

Из формулы (10) следует уравнение динамики механической части:

.

Выразим токи из формул (8) и (9) и получим уравнения динамики электромагнитной части:

;

.

Формулы (20) — (22) и (31) — (33) описывают работу объекта (системы «двигатель? генератор»), которому соответствует структурная схема, приведенная на рисунке 5.

По формулам (4) — (7) с учетом (29) найдем параметры, входящие в уравнения (32), (33):

;

;

;

.

Механический момент внешних сил Mв является суммой нескольких моментов: внешней нагрузки Mвн. н, сухого трения Mтр, сопротивления Mс, зависящего от скорости вращения вала.

Рисунок 5 — Структурная схема модели объекта при включении тяговых двигателей методом взаимной нагрузки В курсовой работе будем считать, что момент внешней нагрузки равен нулю (Mвн.н = 0). Значение момента сухого трения Mтр будем считать равным 0,2 Mн. Действие сухого трения приводит к тому, что вращение двигателя начинается лишь после того, как значение момента вращения вала (Mд? Mг) превысит значение Mтр, поэтому при моделировании следует считать, что при малых значениях (Mд? Mг) возмущение Mтр должно компенсировать полезный момент, а после превышения порогового значения возмущение Mтр остается фиксированной величиной, равной 0,2 Mн. Такой алгоритм можно описать выражением:

.

Динамическая нагрузка, зависящая от скорости вращения вала, описывается формулой:

где в — коэффициент, в = 0,004 Mн.

Таким образом, возмущение Mв состоит из двух составляющих — моментов Mтр и Mс. Для моделирования Mв можно применить подсистему, собранную в пакете Simulink, структура которой представлена на рисунке 6, где переключатель Switch имеет порог переключения, равный малому положительному числу (например, 0,001Щн) (этот блок включает верхний входной сигнал, если скорость больше нуля), а блок насыщения Saturation имеет пороговое значение 0,2Mн.

Рисунок 6 — Модель источника возмущений, действующего на момент вращения вала В соответствии со структурной схемой, представленной на рисунке 5, построим модель объекта в Simulink:

Рисунок 7 — Модель объекта при включении тяговых двигателейметодом взаимной нагрузки Здесь структура блока Subsystem соответствует модели, приведенной на рисунке 6, а блок Fcn реализует нелинейность вида (15).

1.3 Моделирование при программном управлении

Рассмотрим функционирование объекта при одном из типовых входных воздействий. Тяговые двигатели имеют большое значение пускового момента, и при подаче скачка напряжения на их обмотки происходит резкое увеличение тока, многократно превышающего допустимое значение, в связи с этим на практике применяют различные системы запуска, ограничивающие ток и постепенно увеличивающие скорость вращения вала. Таким образом, вместо единичного скачка 1(t) в качестве управляющих напряжений u1(t) и u2(t) будем использовать сигналы, линейно нарастающие до заданного уровня.

Рассмотрим режим работы, когда от нулевого момента времени (t0 = 0) до t1 (3 мин) происходит разгон, затем от t1 до момента t2 (18 мин) поддерживается постоянная скорость, и на отрезке времени от t2 до t3 (20 мин) происходит плавное торможение. Управляющие воздействия на первом отрезке [t0; t1] зададим линейно нарастающими от нуля до значений u1 (1500 В) и u2 (100 В), на втором отрезке [t1; t2]? постоянными (u1 = 1500 В и u2 = 100 В) и на третьем [t2; t3]? линейно убывающими до нуля. Каждое из управляющих воздействий зададим в пакете Simulink с помощью блока Signal Builder, расположенного в библиотеке элементов Sources.

Временные диаграммы изменения управляющих воздействий u1(t) и u2(t) и скорости Щ (t) показаны на рисунках 8 — 10. Временные диаграммы токов iд (t) и iг (t) представлены на рисунке 11.

Рисунок 8 — Временная диаграмма управляющего воздействия u1(t)

Рисунок 9 — Временная диаграмма управляющего воздействия u2(t)

Рисунок 10 — Временная диаграмма изменения скорости Щ (t)

Рисунок 11 — Временная диаграмма изменения токов iд (t) и iг (t)

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

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

Рисунок 12 — Общий вид структурной схемы системы с комбинированным управлением

1.4 Получение модели «вход?сосотояние?выход»

Переменными состояния динамической системы с выходом y (t) называют такие независимые переменные xi (t), которые в каждый момент времени однозначно определяют сигнал y (t).Формой записи дифференциальных уравнений для xi (t) является нормальная (форма Коши), когда уравнения решены или записаны относительно первых производных dxi/dt, а правые части не содержат производных ни от переменных состояния, ни от входного воздействия g (t). В общем случае для многомерных нелинейных динамических систем и объектов нормальная форма Коши записывается виде:

xя (t)=f[x (t), g (t)], (38)

где g (t) — r — мерный вектор входного сигнала, x (t) — вектор переменных состояния размерности n Ч 1.

Выходной сигнал y (t) (вектор размерности k Ч 1) связан с переменными состояния х (t) в общем случае нелинейным соотношением:

y (t)=h[x (t)], (39)

которое так же не содержит производных. Формулу (39)называют уравнением связи с выходом или уравнением наблюдений, так как чаще всего y (t) является измеряемым (наблюдаемым) выходным сигналом объекта или системы. Выражения (38) и (39) называют моделью многомерной динамической системы «вход-состояние-выход». В нашем случае автоматизируется объект управления, поэтому входным сигналом g (t) является управляющее воздействие u (t), описываемое вектором, т. е. r=2,а переменные состояния являются токи двигателя и генератора и скорость.

(40)

(41)

(42)

Перейдем во временную область заменим аргумент s оператором дифференцирования p=d/dt, тогда дифференциальное уравнение для переменной ?(t) будет иметь вид

(43)

Или с учетом выражений (4) и (7)

(44)

Токи в якорных цепях двигателя iд и генератора iг, а также угловая скорость вращения их валов? полностью определяют поведение или динамику электромеханического объекта, поэтому их можно выбрать в качестве переменных состояния:

x1 =?;(45)

x2 = iд;(46)

x3 = iг,(47)

Вектор переменных состояний будет иметь размерность 3Ч1, т. е. n=3. В выбранных координатах выражения (44) запишем уравнение для скорости в виде:

(48)

Те же самые математические операции проведем для токов

(49)

(50)

(51)

(52)

После подстановки введенных обозначений и с учетом формулы (8) получается уравнение для тока двигателя

(53)

То же самое для тока генератора с учетом формулы (9)

(54)

(55)

(56)

(57)

Целью управления является поддержание заданной угловой скоростивращения валов Д и Г, поэтому выходной переменной y (t) является скалярная величина ?, т. е. k=1. В этом случае уравнение связи с выходом (39) становится линейным и записывается следующим образом:, где .

Анализ правых частей дифференциальных уравнений для показывает, что в него входят линейные и нелинейные составляющие от переменных состояния x, а также линейные от управления u и момента Mв, причем последний может рассматриваться как возмущение w (t). Поэтому модель «вход ?состояние?выход» будет иметь вид:

; (58)

(59)

где w — вектор возмущений размерности lЧ1.

В общем случае матрицы A, C, G и имеют соответственно размерности nЧn, nЧr, nЧl и kЧn. Для рассматриваемого объекта управления n=3, r=2,l=1,k=1.

Уравнение состояния имеет вид:

(60)

(61)

1.5 Модели устройства управления

Электромеханический объект управляется двумя напряжениями u 1(t) и u2(t). В случае комбинорованного управления одно из управлений u1(t) заменяется программной траекторией. Временные диаграммы изменения управляющего воздействия u1(t) и скорости ?(t) показаны на рисунке 13.

Рисунок 13 — Временные диаграммы управляющего воздействия и скорости вращения вала при программном управлении При управлении реальным объектом могут возникать отклонения реальных значений параметров объекта от параметров его модели, поэтому при программном управлении невозможно получить точное желаемое значение выходной координаты. Таким образом, необходимо применять принцип управления по отклонению, когда задается желаемое изменение скорости вращения g (t) и определяются управляющие воздействия u1(t) и u2(t) так, чтобы в каждый момент времени минимизировать ошибку:

е (t)=g (t)??(t).(62)

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

Рисунок 14 — Общий вид структурной схемы системы с комбинированным управлением В этом случае модель переменных состояния (58) можно записать в виде:

(63)

(64)

(65)

С учетом перехода к новым обозначениям модель «вход-состояние-выход» имеет вид:

(66)

Здесь g1 — скалярная заданная или программная траектория, а управление u=u2, поэтому матрица С устанавливает связь переменных состояния с u2(t), а G1 c воздействием u1=g1. Размерности матриц С и G1 равны 3Ч1. Принято обозначение G1, так как программная траектория g1 является детерминированным возмущением для моделей «вход-состояние-выход».

В качестве устройства управления используется одномерный ПИД-регулятор с передаточной функцией вида:

(67)

Схема имитационного моделирования ПИД-регулятора в Simulink приведена на рисунке 15.

Рисунок 15 — Схема моделирования ПИД-регулятора Передаточную функцию Wp (s) можно представить следующим образом:

Wp (s)=kp+W (s). (68)

Для выражения W (s) справедливо:

(69)

Степень полинома R (s) числителя m=2 и полинома Q (s) знаменателя n=2, т. е. m=n. В этом случае рекомендуется в функции W (s) выделять целую часть kц и записывать в виде:

(70)

где не известными являются параметр kц и полином R1(s) .

Очевидны формулы:

(71)-(72)

Подставим в последнее равенство выражения для полиномов R (s), Q (s), приведем подобные слагаемые и получим:

(73)

Из выражения (73) следует, что полином R1(s) будет иметь степень m=1, если выполняются равенства:

(74)

тогда (75)

В этом случае передаточная функция ПИД-регулятора записывается следующим образом:

(76)

где л=kp+kd/ф Передаточная функцияW1(s) будет иметь вид:

(77)

где ф1=kiф-kd/ф.

Модели W1(s) соответствует неоднородное дифференциальное уравнение

(78)

Представим выражение (78) в нормальной форме Коши для вектора переменных состояния v, записанного через фазовые координаты u (t) и uя (t):

(79)

Тогда для ПИД-регулятора модель вход-состояние-выход будет иметь вид:

(80)

(81)

где А1- системная матрица размером 2Ч2, С1 и матрица связи с входом е и выходом u размерностей 2Ч1 и 1Ч2 соответственно.

Примем следующие условные обозначения

(82)

(83)

(84)

(85)

Тогда

(86)

(87)

(89)

(90)

(91)

(92)

Модель «вход-состояние-выход» ПИД-регулятора:

(93)

1.6 Методика определения параметров регулятора и модели замкнутой системы

Метод настройки регулятора по ограничению по ошибке основан на оптимизации переходной характеристики h (t) системы, когда параметры регулятора подбираются такими, чтобы переходная характеристика вошла в заранее заданные границы. Для системы «двигатель-генератор» вместо переходной характеристики применяют линейно нарастающие сигналы, поэтому оптимизацию h (t) заменим минимизацией динамической ошибки. Для минимизации следует использовать блок Signal Constraint (NCD Blockset), на вход которого подается сигнал ошибки. Ограничения, накладываемые на сигнал в окне блока Signal Constraint, задаются в диапазоне от -0,01Щн до 0,01Щн, где Щн — номинальная скорость вращения вала.

Зададим в программе m-файла начальные значения коэффициентов:

Ki=-500(94)

Kp=-100(95)

Kd=-10(96)

В качестве настраиваемых параметров нужно добавить коэффициенты kп, kи, kд.

Схема имитационного моделирования в Simulink приведена на рисунке 16.

Рисунок 16 — Схема имитационного моделирования замкнутой системы с ПИД-регулятором Управляющие воздействия u1(t) и u2(t) задаются с помощью блока Signal Builder, входящего в библиотеку элементов Sources. Блок Subsystem моделирует моменты, действующие на валы электромеханической системы.

Рисунок 17 — Настройка регулятора с помощью блока Signal Constraint

В результате настройки регулятора были найдены следующие коэффициенты:

Рисунок 18 — Результат настройки регулятора На рисунках 19 и 20 показаны графики тока двигателя и тока генератора.

Рисунок 19 — Временная диаграмма тока iг

Рисунок 20 — Временная диаграмма тока iд

На рисунке 21 показан график номинальной скорости на выходе системы.

Рисунок 21 — Временная диаграмма скорости вращения вала Выходная переменная y (t) определяется выражением (38), поэтому вместо формулы (62) можно записать:

(97)

В соответствии с структурной схемой (рисунок 16) сигнал е (t) подается на вход ПИД-регулятора, для которого справедливы модель (80) относительно переменных состояния v (t) и уравнение (81) связи с выходным сигналом управления u (t). После подстановки в выражения (80) и (82) формулу (97) вместо е получим для ПИД-регулятора следующие уравнения:

(98)

(99)

Управляющее воздействие u (t)подается на вход объекта, описываемого моделью (63) для переменных состояния x (t).После подстановки в уравнение (63)формулы (99) и приведения подобных, получим:

(100)

(101)

(102)

(103)

(104)

(105)

(106)

и получим на основе выражений (100), (98) и (38) для переменных состояния x (t), v (t) и выходного сигнала y (t) модель вида:

(107)

(108)

(109)

Координаты x и v объединяют одним вектором состояния

(110)

Напомним, что элементы соответствуют следующим физическим величинам:

(111)

Тогда для нелинейной многомерной замкнутой системы управления динамическая модель «вход?состояние ?выход» запишется следующим образом

(112)

(113)

где системная матрица, А имеет размер 5Ч5. В данном случае матрицы связи B, G, G1 со входом g, программной траектории g1 и возмущением w соответственно имеют размер 5Ч1. Матрица связи с выходом

(114)

Так как наблюдается только одна выходная координата z1.

В итоге модель «вход-состояние-выход» объекта и ПИД-регулятора выглядит следующим образом:

(115)

(116)

1.7 Методы численного интегрирования обыкновенных дифференциальных уравнений

Выражение (18) является основной формой представления обыкновенных дифференциальных уравнений. По аналогии запишем модель системы (112), но относительно n-мерной векторной переменной y (t) с элементами yi (t) при i =1,n, учитывая, что выходным сигналом системы является угловая скорость вращения валов двигателей ?(t) и y (t) = ?(t) = y1(t). Переменные g, g1 и w представляют собой заданные функции времени (момент возмущения Mв =0), поэтому замкнутая система управления описывается векторным дифференциальным уравнением:

(117)

С начальными условиями:

(118)

Здесь известными являются начальными условия у0, а неизвестной векторная функция с элементами y1(t), y2(t),…yn (t). В правой части уравнения (117) используется более компактная запись для неавтономных систем, так как g=g (t), g1=g1(t), w=w (t).

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

(119)

Здесь приняты следующие обозначения yk =y (tk), yk+1 =y (tk+1), а под разностью

(120)

понимают шаг интегрирования Подставим разность (119) в уравнение (117) при t=tk и получим:

(121)

Из последнего выражения следует, что для момента времени t = tk+1

численное значение yk+1 может быть определено по формуле:

(122)

которое называют явным методом Эйлера.

Для k =0 при известному у0, вычисляют y1, затем при k=1по полученному y1 определяют y2 и т. д., т. е. алгоритм является рекуррентным.

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

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

(123)

используются значения в середине шага.

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

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

где Tmin? минимальная постоянная времени динамической системы.

При численном интегрировании дифференциальных уравнений (112) постоянная Тmin определяется из выражения

Tmin =min (T1,T2,ф), (125)

где Т1, Т2 — постоянные времени инерциальных звеньев в якорных цепях генератора и двигателя соответственно; ф — постоянная времени звена реального дифференцирования в ПИД-регуляторе.

Составим процедуру вычисления правых частей системы дифференциальных уравнений на Matlab.

function dz=fun1(t, z)

Un=1500; %V

Pn=650 000; %Vt

wn=770*2*pi/60; %ob/min;

KPD=0.927;

Ry=0.0317;

Rv=0.0370;

p=6;

In=Pn/Un/KPD;

Imax=1.2*In;

cE=(Un-In*(Ry+Ry))/In/wn;

cM=Pn/(wn*In²);

L=2*Un/(5*p*wn*In);

J=(6*L*Pn²)/(Ry²*wn²*In²);

Mn=Pn/wn;

beta=0.004*Mn;

R1=2*Rv+Ry;

R2=Ry;

T1=3*L/(2*Rv+Ry);

T2=L/Ry;

K=In;

Kd=-80.57;

Kp=-247.35;

Ki=-474.4;

z1=z (1);

z2=z (2);

z3=z (3);

z4=z (4);

z5=z (5);

if (t<120)

g=75/120*t;

else

if (t>120)&&(t<1080)

g=75;

else

g=-75/120*t;

end

end

if (t<120)

g1=1500/120*t;

else

if (t>120)&&(t<1080)

g1=1500;

else

g1=-1500/120*t;

end

end

Mc=beta*z1;

Md=½tanh (2*z2/Imax)*z2;

Mg=½tanh (2*z2/Imax)*z3;

if ((Md-Mg)<0.2*Mn)&&(z1==0)

Mtr=Md-Mg;

else

Mtr=0.2*Mn;

end

W=-Mtr-Mc;

tau=0.001;

tau1=Ki*tau-Kd*tau;

lambda=Kp+Kd/tau;

funth = tanh (2*z2/Imax);

dz1=-z1/J+(z2-z3)*cM*z2*funth/J-W/J;

dz2=-z2/T1-z1*funth/R1/T1+g1/R1/T1;

dz3=-lambda*z1/R2/T2-z3/T2+z1*funth/R2/T2+lambda*g/R2/T2-g1/R2/T2;

dz4=tau1*z4/tau+z5+tau1*g/tau;

dz5=(Ki/tau-tau1/tau²)*z1-z5/tau+Ki*z5/tau;

dz = [dz1; dz2; dz3; dz4; dz5];

t

[T, z]=ode45(@fun1,[0 1200],[0 0 0 0 0]);

plot (T, z (, 1))

Получим решение модели в переменных состояния стандартной процедурой ode45 при заданных параметрах и начальных условиях.

Рисунок 22 — Полученный график скорости

2. Методы и алгоритмы статистического анализа

Одной из типовых многовариантных процедур в САПР является статистический анализ, к которому относятся методы наихудшего случая и статистических испытаний. В первом случае известны математическая модель проектируемого объекта (замкнутая система управления) в виде уравнений (112), (113) и предельно допустимые отклонения параметров и входных воздействий, при котором выходная переменная y (t) имеет наибольшее отклонение от номинального уровня.

В методе Монте-Карло или статистических испытаний известными являются модель и статистические характеристики в общем случае входных воздействий, параметров и возмущений. В работе предполагается, что случайная величина x (t)?это момент возмущения Mвозм (t), т. е. x (t)=Mвозм (t). Схема моделирования моментов, действующих в электромеханической системе, построена на основе рисунке 4 и с учетом случайной величины Mвозм (t) приведена на рисунке 23.

Рисунок 23 — Схема моделирования моментов Моделью проектируемого объекта является приведенная на рис. 8 схема имитационного моделирования замкнутой системы с ПИД-регулятором.

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

Величина x называется случайной, если в результате испытаний она принимает заранее непредсказуемое значение. Основной математической мо-делью является дифференциальный p (x) или интегральный P (x) законы, которые связаны соотношениями:

(126)

Плотность вероятности p (x) устанавливает связь между значением случайной величины и вероятностью появления этого значения. Функция P (x) характеризует вероятность того, что случайная величина не превзойдет значениях. Нормальный закон распределения определяется следующей плотностью вероятности:

(127)

где m и упараметры закона.

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

(128)-(129)

Моменты связаны с параметрами закона распределения, так в случае нормального закона (127) его параметры m=m1 и у2 =M2.

При моделировании необходимо искусственно получать случайные величины, распределенные по заданному закону распределения. В разделе Sources библиотеки Simulink случайные сигналы генерируются блоком Uniform Random Number, является входным для случайного сигнала Mв (t), т. е. момента возмущения.

2.2 Статистический анализ случайных величин

После имитационного моделирования проводится статистическая обработка полученных результатов, т. е. сигналов x (t), y (t) и при этом их можно рассматривать как случайные величины или случайные процессы.

В первом случае статистический анализ сводится к вычислению оценок моментных характеристики определению эмпирических законов распределения или гистограмм. При этом учитывается дискретный характер выборки, например, для xi при i =1,N, где под N понимается объем выборки.

Оценки математического ожидания и дисперсии вычисляются по следующим формулам:

(130)

Напомним, что величину у называют среднеквадратическим отклонением.

В Simulink для приведения статистического анализа случайных величин используется блок To Workspace, который осуществляет регистрацию сигналов и передает информацию в рабочую область. Схема моделирования приведена на рисунке 24.

Рисунок 24 — Схема моделирования и статистического анализа случайных величин Рисунок 25 — Параметры блока Uniform Random Number (возмущение) Моделирование проведено до момента времени 1200 сек. График случайного сигнала показан на рисунке 26:

Рисунок 26 — График входного возмущающего воздействия Статистические данные, полученные в Маtlab:

datastats (Mvoz.signals.values)

ans =

num: 1201

max: 0.9989

min: -0.9999

mean: -0.0087

median: -0.0068

range: 1.9988

std: 0.5653

На рисунке 27 представленная полученная эмпирическая гистограмма для сигнала возмущения, которая соответствует нормальному закону распределения.

Рисунок 27 — Эмпирическая гистограмма для сигнала Мвозм

2.3 Статистический анализ случайных процессов

Если случайные возмущения x (t) =Mвозм и выходную переменную y (t) динамической системы рассматривать как функции времени, то применяют методы статистического анализа случайных процессов. Функция или процесс называется случайным, если ее мгновенные значения в любые дискретные моменты времени являются случайными величинами. Полной характеристикой случайного процесса x (t) является закон распределения p (x, t). Процессы, у которых плотности распределения инвариантны относительно t, т. е. p (x, t) = p (x), называются стационарными. Для стационарных эргодических процессов математическое ожидание и дисперсия являются постоянными величинами. Статистический анализ случайных процессов проводят во временной областии на основе корреляционных функций, а в комплексной (частотной) — на основе спектральных характеристик. Во временной области вводится понятие смешанной моментной характеристики второго порядка. Она называется автокорреляционной функцией и определяет статистическую связь между мгновенными значениями процесса x (t) в разные моменты времени. Для стационарных эргодических процессов она имеет вид:

(131)

Взаимные корреляционные функции оценивают тесноту связи между случайными процессами x (t) и y (t) и определяются по формуле:

(132)

В формулах (131) и (132) под m1(x) и m1(y) обозначены математические ожидания процессов x (t) и y (t) соответственно. По корреляционным функциям можно определить максимальное время корреляции фмк выполняется. Таким образом, если при дискритизации случайных процессов выбрать шаг, то дискретная выборка должна состоять из независимых (некоррелированных) значений.

При проведении статистического анализа определяют выборочные оценки корреляционных функций, которые получают из соответствующих формул (101) или (102) при замене операций интегрирования на суммирование. Например, для случайного процесса оценка имеет вид:

(133)

где ?t — период дискретизации; N — объем выборки; m — оценка среднего значения процесса x (t).

В состав пакета Simulink входит блок Cross-Correlator, с помощью которого можно вычислять как автокорреляционные, так и взаимные корреляционные функции. Схема имитационного моделирования системы и статистического анализа во временной области приведена на рисунке 28.

Рисунок 28 — Схема имитационного моделирования системы и статистического анализа во временной области управление статистический система моделирование В частотной области для случайных процессов определены спектральные плотности мощности, или сокращенно спектральные плотности, которые связаны с корреляционными функциями Rx (ф) и Rxy (ф) следующими выражением:

(134)-(135)

т.е. Sx (щ) и Sxy (щ) являются прямыми преобразованиями Фурье от Rx (ф) и Rxy (ф) соответственно.

Справедливы соотношения:

(136)

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

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

(137)

где ?щ — выбранный шаг дискретизации по частоте.

Формула (137) содержит оценки корреляционной функции и поэтому точность вычислений спектральных плотностей часто недостаточно высока. Поэтому в состав Simulink помимо блока Power Spectral Density, предназначенного для определения оценок спектральной плотности мощности входит и блок Avarging Power Spectral Density, осуществляющей и усреднение оценок.

Преобразование Фурье можно применить и непосредственно к случайному процессу x (t):

(138)

и в этом случае говорят о спектре сигнала.

Вычисление оценок x (щ) осуществляется блоками и анализаторами спектра Spectrum Analyses и Avarging Spectrum Analyses, в последнем производится и усреднение оценок.

Рисунок 29 — Автокорреляционная функция возмущения Rx (ф) По рисунку 29 видно, что сигнал x (t) является некоррелированным.

Рисунок 30 — Спектральная плотность Sx (щ) По рисунку 30 видно, что сигнал x (t) является белым шумом, так как спектральная плотность равномерно распределена на частотном интервале от 0 до 30 рад/сек.

Рисунок 31 — Автокорреляция скорости Ry (ф) Рисунок 32 — Спектральная плотность скорости Sy (щ) На рисунке 32 видно, что спектральная плотность сосредоточена в области низких частот.

Рисунок 33 — Взаимнокорелляционная функция Rxy (ф) В результате получены автои взаимно корреляционные функции, по которым видно что сигнал x (t) является некоррелированным.

Заключение

В ходе курсового проектирования построена модель объекта автоматизации — диагностического стенда испытаний тяговых двигателей. Проведено моделирование системы при программном управлении, когда на вход подаются программно заданные управления u1 и u2. В результате моделирования получены временные диаграммы скорости вращения вала и токов в обмотках двигателя и генератора. Установившиеся значения скорости и токов превысили номинальные значения, в связи с чем было принято решение о введении в систему регулятора.

По передаточной функции регулятора, используя аппарат переменных состояний, была получена модель «вход? состояние? выход». Также составлена схема имитационного моделирования ПИД-регулятора в Simulink. Также была построена в Simulink схема имитационного моделирования замкнутой системы и осуществлена настройка (определение численных значений параметров) регулятора по минимуму динамической ошибки в системе. Была промоделирована работа системы и построены временные зависимости для управляющих напряжений u1 и u2 и выходной переменной y (t)? скорости вращения (t).

Были записаны математические модели ПИД-регулятора и объекта относительно входного сигнала g (t) и вектора состояния x (t). Получена модель «вход? состояние? выход» для расширенного вектора состояния z (t) = |x (t) v (t)|? замкнутой системы управления. Проведено имитационное моделирование и статический анализ случайных величин. Были получены статистические характеристики для сигналов x (t) и y (t). В результате сравнения численных значений математического ожидания и дисперсии, полученных с помощью команд dataststs и var сделан вывод, что они согласуются с заданными параметрами блока Random Number. Для данного числа интервалов группирования была получена и построена эмпирическая гистограмма, соответствующая нормальному закону распределения. Проведено имитационное моделирование и статистический анализ случайных процессов во временной области и спектральных характеристик случайных процессов, также получены и построены оценки автокорреляционных функций Rx (ф) и Ry (ф). В результате сравнения полученных характеристик с идеальной характеристикой белого шума, был сделан вывод, что сигнал x (t) является белым шумом, так как спектральная плотность равномерно распределена на частотном интервале от 0 до 30 рад/сек. В результате получены автои взаимнокорреляционные функции, по которым видно что сигнал x (t) является некоррелированным (рисунки 29−33).

Библиографический список

1. Лаврухин А. А., Когут А. Т. Автоматизированное проектирование нелинейных систем управления: методические указания к курсовому проектированию по дисциплине «Автоматизация проектирования систем и средств управления»; ОмГУПС. Омск, 2011, 34 с.

2. Бесекерский В. А., Попов Е. П. Теория систем автоматического управления, СПб, 2003, 752 с.

3. СТП ОмГУПС -1.2−2005 Работы студенческие учебные и выпускные квалификационные. Общие требования и правила оформления текстовых документов: стандарт предприятия. Омск, 2005, 27 с.