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

Оглавление

• Часть 1. Анализ объекта управления
• 1. Постановка задачи
• 2. Математические модели объекта управления
• 2.1 Уравнение в переменных состояния
• 2.2 Передаточная функция
• 2.3 Весовая функция
• 2.4 Уравнение вход-выход
• 2.5 Частотные характеристики
• 3. Свойства системы
• 3.1 Устойчивость
• 3.2 Анализ минимально фазовости объекта
• 3.3 Исследование управляемости и наблюдаемости
• 3.4 Анализ установившихся режимов
• 3.5 Окончательный выбор параметров и его обоснование
• 4. Процессы в объекте управления
• 4.1 Импульсное воздействие
• 4.2 Ступенчатое воздействие
• 4.3 Гармоническое воздействие
• Часть 2. Синтез законов управления для систем с обратной связью
• 1. Структурная схема системы с регулятором
• 2. Настройка контура управления
• 3. Настройка контура оценивания
• 4. Завершение построения системы
• 5. Сравнение результатов автоматического управления по средством обратнай связи с командным управлением
• Выводы
• Приложение

Часть 1. Анализ объекта управления

1. Постановка задачи

В данной работе рассматривается модель развития многоотраслевой экономики В.В. Леонтьева в частном случае для двух отраслей.

Выбранные данные для модели:

2. Математические модели объекта управления

2.1 Уравнение в переменных состояния

Вектор состояния представим в следующем виде:

Тогда система линейных дифференциальных уравнений, описывающих систему, выглядит

соответственно запишем следующие матрицы и векторы:

2.2 Передаточная функция

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

Передаточная функция, вычисленная при помощи символьной алгебры в MatLab по той же формуле, которая совпадает с найденной аналитически (М-файл №1 в приложении):

H = (k* (10*k1 + a2*k2 + k1*p)) / (p* (p^2 + 110*p - a1*a2 + 1000)) - ( (k - 1)

* (100*k2 + a1*k1 + k2*p)) / (p* (p^2 + 110*p - a1*a2 + 1000))

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

Очевидно, выполняется равенство, которое следует из блок-схемы и структурных свойств систем управления:

2.3 Весовая функция

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

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

, ,

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

Для которой, обратное преобразование Лапласа можно провести с помощью таблицы:

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

Решая систему, получаем соответственно следующие значения коэффициентов:

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

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

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

h =( (cosh (t* (a1*a2 + 2025) ^ (1/2)) + (sinh (t* (a1*a2 + 2025) ^ (1/2)) * ( (10000*k2 + 110*a1*k1 + 100*k*k1 - 10000*k*k2 + a1*a2*k2 - 110*a1*k*k1 + 110*a2*k*k2 + a1*a2*k*k1 - a1*a2*k*k2) / (100*k2 + a1*k1 + 10*k*k1 - 100*k*k2 - a1*k*k1 + a2*k*k2) - 55)) / (a1*a2 + 2025) ^ (1/2)) * (100*k2 + a1*k1 + 10*k*k1 - 100*k*k2 - a1*k*k1 + a2*k*k2)) / (exp (55*t) * (a1*a2 - 1000)) - (100*k2 + a1*k1 + 10*k*k1 - 100*k*k2 - a1*k*k1 + a2*k*k2) / (a1*a2 - 1000)

Сравним теоретические результаты с результатами полученными численно, с помощью MatLab, подставив определенные численные значения и сравним их (М-файл №2 в приложении):

h1 = 0.0217

h = 0.0217

hm= 0.0217

Где h1 - ответ, полученный аналитически

h - ответ, полученный при обратном преобразовании Лапласа

hm - ответ, полученный при вычислении через переменные состояния

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

2.4 Уравнение вход-выход

Обозначим числитель и знаменатель передаточной функции согласно формуле:

Тогда уравнение вход-выход запишется в следующем виде, если вместо переменной р подставить оператор дифференцирования по времени:

2.5 Частотные характеристики

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

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

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

Для символьного решения посредством MatLab ограничимся вычислением частотных характеристик через переменные состояния, чтобы операции не сводились к простому переименованию переменных. Проведем проверку между аналитическим и символьным расчетом в MatLab, при фиксированных значениях параметров (М-файл №3 в приложении):

H = - (k* (10*k1 + a2*k2 + k1*w*i)) / (w* (w^2*i + 110*w + a1*a2*i - 1000*i)) + ( (k - 1) * (100*k2 + a1*k1 + k2*w*i)) / (w* (w^2*i + 110*w + a1*a2*i - 1000*i))

modh = 0.0107

modh1 =0.0107

argh =-1.7168

argh1 =-1.7168

Ответы совпали, значит, аналитические вычисления верны.

3. Свойства системы

3.1 Устойчивость

Найдем собственные числа матрицы А (очевидно это корни знаменателя передаточной функции)

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

Для того, чтобы удовлетворить условию Стодолы, которое является в нашем случае не только необходимым условием, но и достаточным, потребуем следующее:

,

Корневой годограф смотри в приложении.

3.2 Анализ минимально фазовости объекта

Запишем уже известные нам передаточные функции:

,

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

очевидно, неравенство (3) выполняется, если выполнены первые два, тогда остаются

3.3 Исследование управляемости и наблюдаемости

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

Корни знаменателя мы знаем из исследования устойчивости:

,

корни числителей:

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

При этом стоит ввести вместо строгого отсутствия равенства определенное неравенство, чтобы избежать слабоуправляемых и плохо наблюдаемых систем:

Соответственно запишем, получившиеся ограничения на наши коэффициенты:

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

Теперь будем исследовать управляемость и наблюдаемость посредством MatLab. Составим матрицу управления и матрицу наблюдаемости по формулам через переменные состояния:

Полученные решения:

P =

[k, - 100*k - a1* (k - 1), k* (a1*a2 + 10000) + 110*a1* (k - 1)]

[1 - k, 10*k + a2*k - 10, - 110*a2*k - (a1*a2 + 100) * (k - 1)]

[0, k*k1 - k2* (k - 1), (10*k2 - a1*k1) * (k - 1) - k* (100*k1 - a2*k2)]

Q =

[0, conj (k1), conj (a2) *conj (k2) - 100*conj (k1)]

[0, conj (k2), conj (a1) *conj (k1) - 10*conj (k2)]

[1, 0, 0]

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

Находим в MatLab символьные выражения для определителей матриц управления и наблюдаемости (М-файл №4 в приложении):

dP = 9000*k-100*a1+380*a1*k+190*k^2*a2-17100*k^2-460*a1*k^2-180*k^3*a2+180*a1*k^3-k*a1*a2+3*a1^2*k-3*a1^2*k^2+a1^2*k^3+k^3*a2^2-a1^2+8100*k^3+3*k^2*a1*a2-2*k^3*a2*a1

dQ = conj (a1) *conj (k1) ^2 + 90*conj (k1) *conj (k2) - conj (a2) *conj (k2) ^2

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

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

При найденных значениях (см. раздел 3.5) параметров:

dP = 518.4000

dQ = 99.1500

Что значительно больше любого , которое на один-два порядка меньше левой части.

3.4 Анализ установившихся режимов

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

3.5 Окончательный выбор параметров и его обоснование

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

1)

2)

3)

4)

5)

6)

7)

8)

9)

10)

11)

12)

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

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

Перепишем первое уравнение:

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

Очевидно также то, что при выполнении условия (5) (физический смысл параметров) автоматически выполняются условия (2) и (3). Таким образом, система неравенств для выглядит следующим образом:

Ограничение (2) позволяет нам и на этом этапе выкинуть некоторые неравенства. Так как есть некоторая небольшая окрестность вокруг точки , в которой левая часть неравенства обращается в ноль, то если эта не является отрицательной, мы можем смело выкинуть это условие. Таким образом, избавляемся от условий (4), (5)

Решая неравенство (1) и объединяя их в систему с неравенством (6), получаем:

Таким образом, мы имеем понятие о порядке левой части неравенств (3) и (7), и можем теперь подобрать , которое будет соответствовать нашим запросам, подставляя значения из промежутка, убеждаемся, что левая часть получается порядка десяток и сотен.

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

Решая систему неравенств получаем итоговый интервал изменения

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

Таким образом, окончательно выбранные параметры системы: , , , ,

4. Процессы в объекте управления

4.1 Импульсное воздействие

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

Строим графики реакции на импульсное воздействие при помощи MatLab (М-файл №5 в приложении):

Сравним данные полученные с помощью MatLab c данными полученными теоритически:

a. Установившийся режим:

b. Начальные точки.

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

4.2 Ступенчатое воздействие

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

Построим реакции, полученные при помощи MatLab (М-файл №5 в приложении):

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

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

И то же самое можно показать и для реакции всей системы

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

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

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

По причине того, что при заданных параметрах мы получаем чисто вещественные корни:

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

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

4.3 Гармоническое воздействие

Для начала опять рассмотрим характеристики инерциальных звеньев. Для этого воспользуемся пакетом MatLab, в котором построим графики их частотных характеристик (М-файл №5 в приложении).

Аналогично могут быть посчитаны и построены эти же графики, но уже аналитическим методом, который предложен в разделе 2.5.

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

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

Часть 2. Синтез законов управления для систем с обратной связью

1. Структурная схема системы с регулятором

В данной работе строим управление системой по оценке вектора состояния. Приведем структурную схему.

И полная схема замкнутой системы:

2. Настройка контура управления

Будем строить управление по оценке вектора состояния:

Где - управление, зависящее от задания желаемого выхода, - оценка вектора состояния. ля начала найдем просто матрицу , которая в нашем случае является строкой. Для этого запишем характеристическое уравнение матрицы A:

,

,

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

Приведем уравнение с к виду:

Таким образом, параметры : и собственные значения матрицы А:

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

,

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

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

Запишем характеристическое уравнение матрицы .

Приведем это уравнение к виду:

При перемножении и подстановки выбранных желаемых собственных чисел:

Вводим промежуточный вектор , который зависит от установленных значений собственных чисел замкнутой системы и исходной по формуле:

А искомый К определяется:

,

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

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

>> k=place (A,b, [-100 - 90 - 80])

1.0e+004 * 0.0152 0.0191 1.9565

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

Построим матрицу :

3. Настройка контура оценивания

Так как , для управления по вектору состояния, нам требуется найти его оценку по известной модели и выходу. Таким образом, строим алгоритм оценивания, в котором требуется матрица невязки L. А - ошибка оценивания.

,

Введем обозначение (отсюда можно заключить, что L вектор):

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

И делаем аналогичные предыдущему разделу операции. Запишем характеристическое уравнение матрицы.

Приведем это уравнение к виду:

При перемножении и подстановки выбранных желаемых собственных чисел:

И уже найденные коэффициенты: .

Вводим промежуточный вектор , которая зависит от установленных значений собственных чисел контура оценивания и исходной матрицы А:

А искомый L определяется:

,

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

, , , ,

,

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

>> l=place (A',c', [-110 - 120 - 130]) l = 1.0e+004 * 0.8940 1.1482 0.0250

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

Построим матрицу :

4. Завершение построения системы

Запишем полученную замкнутую систему в виде:

Для простоты и наглядности получения параметров префильтра перейдем к описанию системы в виде передаточной функции. С легкостью можно получить передаточную функцию системы, подставив в MatLab описание замкнутой системы в переменных состояния (М-файл №6 в приложении):

Transfer function:

0.78s+36.8

-----------------------------------

s^3 + 270 s^2 + 2.42e004 s + 7.2e005

Если с помощью MatLab посчитать корни знаменателя, то они полностью совпадут с желаемыми:

>> roots ([1 270 2.42e004 7.2e005])

100.0000

90.0000

80.0000

5. Сравнение результатов автоматического управления по средством обратнай связи с командным управлением

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

Система без обратных связей:

Система управляемая автоматически:

Выводы

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

Приложение

М-файл №1. Часть 1. Раздел 2.2 (определение передаточной функции).

syms a1 a2 k p k1 k2

A= [-100 a1 0; a2 - 10 0; k1 k2 0];

B= [k; 1-k; 0];

C= [0 0 1];

E= [1 0 0; 0 1 0; 0 0 1];

H=C*inv (p*E-A) *B

М-файл №2. Часть 1. Раздел 2.3 (определение весовой функции).

syms t a1 a2 k p k1 k2

hm=C*expm (A*t) *B

pause

h=ilaplace (k* (10*k1 + a2*k2 + k1*p)) / (p* (p^2 + 110*p - a1*a2 + 1000)) - ( (k - 1) * (100*k2 + a1*k1 + k2*p)) / (p* (p^2 + 110*p - a1*a2 + 1000)),t) (1)

pause

a1=1;

a2=2;

k=0.8;

t=2;

k1=0.8;

k2=0.7;

h=… (2)

hm=… (3)

h1=… (4)

(1) выражение взято из М-файла №1

(2) ответ, полученный в четвертой строчке

(3) ответ, полученный во второй строчке

(4) ответ, полученный аналитически

М-файл №3. Часть 1. Раздел 2.5 (частотные характеристики).

syms a1 a2 k w k1 k2

A= [-100 a1 0; a2 - 10 0; k1 k2 0];

B= [k; 1-k; 0];

C= [0 0 1];

E= [1 0 0; 0 1 0; 0 0 1];

H=C*inv (i*w*E-A) *B

pause

a1=1;

a2=2;

k=0.8;

k1=0.8;

k2=0.7;

w=2;

modh=abs (…) (1)

modh1=… (2)

argh=angle (…) (3)

argh1=… (4)

(1) выражение для Н (шестая строчка)

(2) аналитическое выражение модуля

(3) выражение для Н (шестая строчка)

(4) аналитическое выражение аргумента

М-файл №4. Часть 1. Раздел 3.3 (управляемость и наблюдаемость).

syms a1 a2 k k1 k2

A= [-100 a1 0; a2 - 10 0; k1 k2 0];

B= [k; 1-k; 0];

C= [0 0 1];

E= [1 0 0; 0 1 0; 0 0 1];

P= [B A*B A*A*B]

Q= [C' A'*C' A'*A'*C']

pause

dP=det (P)

dQ=det (Q)

pause

a1=80;

a2=5;

k=0.8;

dP=… (1)

dQ=… (2)

(1) выражение, полученное для dP

(2) выражение, полученное для dQ

М-файл №5. Часть 1. Раздел 4 (импульсное, ступенчатое и гармоническое воздействия).

sys1=tf ([0.8 23], [1 110 410.717]);

sys2=tf ([0.2 26.28568], [1 110 410.717]);

sys3=tf ([0.78 36.8], [1 110 410.717 0]);

impulse (sys1)

grid on

hold on

impulse (sys2)

figure

impulse (sys3)

grid on

figure

step (sys1)

grid on

hold on

step (sys2)

figure

step (sys3)

grid on

figure

bode (sys1)

hold on

grid on

bode (sys2)

figure

grid on

bode (sys3)

figure

nyquist (sys1)

hold on

grid on

nyquist (sys2)

figure

grid on

nyquist (sys3)

М-файл №6. Часть 2. Раздел 4 (Завершение построения системы)

A= [-100 75 0; 7.8571 - 10 0; 0.8 0.7 0];

B= [0.8; 0.2; 0];

C= [0 0 1];

K=10000* [0.0152 0.0191 1.9565];

L=10000* [0.8940 1.1482 0.0250];

A1= [A-B*K B*K; zeros (3,3) A-L'*C];

B1= [B; 0; 0; 0];

C1= [C 0 0 0];

sys=ss (A1,B1,C1,0)

Корневой годограф.

q=-4000: 5: 4000;

for i=1: 1601

l1 (i) =-55-sqrt (2025+q (i));

l2 (i) =-55+sqrt (2025+q (i));

plot (real (l1 (i)), imag (l1 (i)),real (l2 (i)), imag (l2 (i)))

hold on

end

Размещено на Allbest.ru