Цифровые регуляторы

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

Техническое задание

Рис 1. Схема системы

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

1. Построить переходную характеристику объекта, реализовать кривую разгона y(t) и выбрать по ней интервал квантования Т₀.

2. Произвести имитационное моделирование объекта управления с возмущающим воздействием, произвести регистрацию реакции возмущающего воздействия с заданной корреляционной функцией:

R_n(ф)= ;

3. Построить дискретную модель объекта управления следующими методами:

· Переходом от дифференциального уравнения к разностному;

· Переходом от непрерывной передаточной функции W₀(p) объекта к дискретной G₀(z) по таблицам z-преобразования;

· Построением дискретно совпадающей модели объекта по кривой разгона;

· Построением МНК-модели по кривой разгона;

Выбрать адекватную дискретную модель из всех построенных.

4. Произвести имитационное моделирование цифровой управляемой системы с цифровым обобщенным линейным регулятором второго или более высокого порядка. Произвести параметрическую оптимизацию цифрового регулятора методом Гаусса-Зайделя по критерию минимума квадратичной интегральной оценки:

J= = min;

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

5. Рассчитать параметры цифрового регулятора с минимальной обобщенной дисперсией и произвести имитационное моделирование соответствующей управляемой системы.

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

7. Сделать рекомендации по выбору вида цифрового регулятора для заданного объекта управления.

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

W₀(p)= , где с = 0,5; K₀ = 2; T = 0,2 c; D=0,3;

Подставим значения параметров и получим передаточную функцию рассматриваемого объекта управления:

W₀(p)= =;

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

R_n(ф)= ;

=== 5; === 1;

= 2; = 4;

Используя полученные исходные данные выполним техническое задание.

1. Имитационное моделирование объекта управления

С помощью пакета «MATLAB» построим переходную характеристику рассматриваемого объекта управления. Как было рассчитано выше, передаточная функция объекта управления имеет следующий вид:

W₀(p)= ;

Для получения переходной характеристики объекта управления введем следующую команду «MATLAB»:

>> step(tf([1],[0.01, 0.06, 1]));

Получим следующее изображение переходной характеристики:

Рис 1.1. Переходная характеристика объекта управления с квантованием по времени

Согласно данному изображению переходной характеристики выберем интервал квантования Т₀ = 0,1 c.

2. Имитационное моделирование возмущающего воздействия и построение эмпирической оценки его корреляционной функции

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

Рис 3.1. Формирующий фильтр с «белым шумом» на входе

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

Рис 3.2. Параллельное соединение двух формирующих фильтров

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

R_n(ф)= ;

Подставим в данную корреляционную функцию исходные данные:

R_n(ф)= = 2+ 4 = R_n1(ф)+ R_n2(ф)

Построим для наглядности графики корреляционных функций в системе «MATLAB»:

>> syms t;

>> ezplot(2*exp(-5*abs(t)))

>> ezplot(4*exp(-1*abs(t)))

>> ezplot(4*exp(-1*abs(t))+ 2*exp(-5*abs(t)))

Рис 3.3. Корреляционные функции отдельных составляющих блоков формирующих фильтров



Рис 3.4. Корреляционная функция формирующего фильтра

По известной корреляционной функции R_nn(ф) найдем спектральную плотность возмущающего воздействия S_nn(щ), характеризующую частотный состав процесса и определяющую распределение среднего значения мощности по спектру:

;

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

.

Согласно определению белого шума его спектральная плотность является постоянной величиной. Пусть S_v(щ) = const = 1, тогда:

;

;

;

;

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

= 0,9; = 0,2; ;

= 2,8; = 1; ;

W_фф(p) = + ;

Учитывая выбранный ранее шаг квантования (Т₀ = 0,1 с), определим теперь дискретную передаточную функцию формирующего фильтра, используя следующую замену:

;

;

= ;

= ;

G_фф(z) = +;

Построим в «Simulink» дискретную модель формирующего фильтра и снимем временные характеристики белого и окрашенного шума:

Рис 3.5. Модель формирующего фильтра



Рис 3.6. Белый шум

Рис 3.7. Окрашенный шум

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

3.1 Построение дискретной модели переходом от дифференциального уравнения к разностному

Передаточная функция рассматриваемого объекта управления имеет следующий вид:

W₀(p)= ==> (0,01p²+0,06p+1)·Y(p)=U(p);

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

0,01+0,06+y(t) = u(t);

Перейдем теперь от дифференциального уравнения к разностному, используя левую разность и учитывая размер такта Т₀ = 0,1 c:

Д²y(k) = y(k) - 2y(k-1) + y(k-2), Дy(k) =y(k) - y(k-1) =>

+ + y(k) = u(k);

+ + y(k) = u(k);

y(k) - 2y(k-1) + y(k-2) + 0,6y(k) - 0,6y(k-1)+y(k)=u(k) =>

2,6y(k) - 2,6y(k-1) + y(k-2) = u(k) => y(k) = =>

y(k)=y(k-1) - 0,4y(k-2) + 0,4u(k);

В полученном выражении y(k) - выходной сигнал, u(k) - входной сигнал (единичное ступенчатое воздействие), k - отсчеты времени. C помощью данного выражения рассчитаем теперь значения на кривой разгона в моменты квантования, учитывая что график на рисунке 1.1 выходит из начала координат:

y(0)=0;

y(1)=y(Т₀)=y(0) - 0,4y(-1) + 0,4u(1)= 0 - 0 + 0,4 = 0,4;

y(2)=y(2Т₀)=y(1) - 0,4y(0) + 0,4u(2)= 0,4 - 0 + 0,4 = 0,8;

y(3)=y(3Т₀)=y(2) - 0,4y(1) + 0,4u(3)= 0,8 - 0,16 + 0,4 = 1,04;

y(4)=y(4Т₀)=y(3) - 0,4y(2) + 0,4u(4)= 1,04 - 0,32 + 0,4 = 1,12;

y(5)=y(5Т₀)=y(4) - 0,4y(3) + 0,4u(5)= 1,12 - 0,42 + 0,4 = 1,1;

y(6)=y(6Т₀)=y(5) - 0,4y(4) + 0,4u(6)= 1,1 - 0,45 + 0,4 = 1,05;

y(7)=y(7Т₀)=y(6) - 0,4y(5) + 0,4u(7)= 1,05 - 0,44 + 0,4 = 1,01;

y(8)=y(8Т₀)=y(7) - 0,4y(6) + 0,4u(8)= 1,01 - 0,42 + 0,4 = 0,99;

y(9)=y(9Т₀)=y(8) - 0,4y(7) + 0,4u(9)= 0,99 - 0,4 + 0,4 = 0,99;

y(10)=y(10Т₀)=y(9) - 0,4y(8) + 0,4u(10)= 0,99 - 0,4 + 0,4 = 0,99;

Чтобы получить передаточную функцию, применим к выражению для y(k) z-преобразование:

Z{y(k)} = Z{y(k-1) - 0,4y(k-2) + 0,4u(k)};

Y(z) = Y(z)·z^-1 - 0,4·Y(z)·z^-2 +0,4·U(z) =>

Y(z)·(1 - z^-1 + 0,4·z^-2) = 0,4·U(z) =>

G₀(z)= =;

Рассчитаем статический коэффициент передачи объекта, подставив вместо z единицу:

K₀ =G₀(1)= =1;

3.2 Построение дискретной модели переходом к дискретной передаточной функции объекта

Перейдем от непрерывной передаточной функции объекта управления W₀(p) к дискретной G₀(z), используя таблицы z-преобразования:

G₀(z)=(1-z^-1)·Z =(1-z^-1)·Z =(1-z^-1)·Z ;

Для перехода воспользуемся следующей зависимостью, учитывая величину такта квантования Т₀ = 0,1 c:

p^-1= = 0,05

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

G₀(z)= ==

= = ==

= =

==;

Таким образом, в итоге получаем следующее выражение для дискретной передаточной функции G₀(z):

G₀(z) = ;

Рассчитаем теперь статический коэффициент передачи объекта, подставив вместо z единицу:

K₀ =G₀(1)= = = 1;

3.3 Построение дискретно-совпадающей модели по кривой разгона

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

G₀(z)= == =>

Y(z)·( 1 + a₁·z^-1 + a₂·z^-2) = U(z) ·( b₀ + b₁·z^-1);

Y(z)+Y(z)·a₁·z^-1 +Y(z)·a₂·z^-2 = U(z)·b₀ + U(z)·b₁·z^-1;

Применим к данному выражению обратное z-преобразование:

y(k) + a₁·y(k-1)+ a₂·y(k-2) = b₀ ·u(k) + b₁·u(k-1);

y(k) = b₀ ·u(k) + b₁·u(k-1) - a₁·y(k-1) - a₂·y(k-2);

В данном выражении y(k) - выходной сигнал, u(k) - входной сигнал (единичное ступенчатое воздействие), k - отсчеты времени. Из рисунка 1.1 видно, что кривая разгона выходит из начала координат, т.е. y(0)=0. Следовательно:

y(0) = b₀ ·u(0) + b₁·u(-1) - a₁·y(-1) - a₂·y(-2)= b₀ ·1+ 0 - 0 - 0=0 => b₀ =0;

Исходя из полученного условия (b₀ =0), выражение для y(k) можно записать следующим образом:

y(k) = b₁·u(k-1) - a₁·y(k-1) - a₂·y(k-2);

Согласно выбранному интервалу квантования Т₀ = 0,1 c возьмем несколько точек на графике, изображенном на рисунке 1.1:

y(0) = y(0Т₀) = 0;

y(1) = y(Т₀) = 0,39;

y(2) = y(2Т₀) = 1,02;

y(3) = y(3Т₀) = 1,35;

y(4) = y(4Т₀) = 1,29;

y(5) = y(5Т₀) = 1,06;

y(6) = y(6Т₀) = 0,89;

y(7) = y(7Т₀) = 0,87;

y(8) = y(8Т₀) = 0,95;

y(9) = y(9Т₀) =1,03;

y(10) = y(10Т₀) = 1,05;

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

y(3) = b₁·u(2) - a₁·y(2) - a₂·y(1) = b₁ - a₁·1,02 - a₂·0,39 = 1,35;

y(6) = b₁·u(5) - a₁·y(5) - a₂·y(4) = b₁ - a₁·1,06 - a₂·1,29 = 0,89;

y(10) = b₁·u(9) - a₁·y(9) - a₂·y(8) = b₁ - a₁·1,03 - a₂·0,95 = 1,05;

Решив систему данных уравнений, получим следующие значения неизвестных параметров:

a₁ = -0,93; a₂ = 0,55; b₁ = 0,62;

Подставим полученные коэффициенты в исходное выражение для y(k) и получим дискретно-совпадающую модель:

y(k) = 0,62·u(k-1) + 0,93·y(k-1) - 0,55·y(k-2);

Подставив найденные выше коэффициенты еще и в выражение для дискретной передаточной функции, получим следующее выражение:

G₀(z) = ;

Рассчитаем теперь статический коэффициент передачи объекта, подставив вместо z единицу:

K₀ =G₀(1)= = = 1;

3.4 Построение МНК-модели по кривой разгона

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

y(k) = b₁·u(k-1) - a₁·y(k-1) - a₂·y(k-2);

Найдем теперь неизвестные коэффициенты в данном выражении методом наименьших квадратов. Представим данное выражение в матричном виде:

y = V·в;

В данном выражении y - вектор экспериментально снятых отсчетов, V - матрица входов-выходов, в - вектор неизвестных параметров модели:

y = = ; V = = ; в =

;

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

в = (V^T·V) ^-1 · V^T ·y;

Для решения данного уравнения воспользуемся пакетом «MATLAB», введя следующие команды:

>> y=[0.39; 1.02; 1.35; 1.29; 1.06; 0.89; 0.87; 0.95; 1.03; 1.05]

>> V=[1 0 0; 1 0.39 0; 1 1.02 0.39; 1 1.35 1.02; 1 1.29 1.35; …

… 1 1.06 1.29; 1 0.89 1.06; 1 0.87 0.89; 1 0.96 0.87; 1 1.03 0.95];

>> B=((V'*V)^(-1))*V'*y;

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

в = ;

Тогда неизвестные в выражении для y(k) имеют следующие значения:

b₁ = 0,4934; a₁ = -1,086; a₂ = 0,5954;

Подставив данные значения в исходное уравнение, получим выражение для МНК-модели:

y(k) = 0,4934·u(k-1) - 1,086·y(k-1) + 0,5954·y(k-2);

Подставив найденные выше коэффициенты еще и в выражение для дискретной передаточной функции, получим следующее выражение:

G₀(z) = ;

Рассчитаем теперь статический коэффициент передачи объекта, подставив вместо z единицу:

K₀ = G₀(1)= = 0,97;

3.5 Сравнение полученных моделей

Построим переходные процессы полученных дискретных моделей и сравним их с исходным переходным процессом. Для этого используем следующие команды «MATLAB»:

>> W = tf([1],[0.01, 0.06, 1]);

>> G1 = tf([0.4],[1 -1 0.4], 0.1);

>> G2 = tf([1.25 3.75 3.75 1.25],[15.5 -15 9.5], 0.1);

>> G3 = tf([0.62 0],[1 -0.93 0.55], 0.1);

>> G4 = tf([0.4934 0],[1 -1.086 0.5954], 0.1);

>> step(W, G1, G2, G3, G4);



Рис 2.5.1. Переходные процессы моделей

Из рис. 2.5.1 очевидно, что наиболее адекватной является четвертая модель, т.е. МНК-модель системы управления. В дальнейшем будем использовать ее.

4. Настройка системы с цировым ПИД-регулятором

4.1 Настройка ПИД-регулятора непрерывной модели системы

Рис 4.1.1. Схема непрерывной модели системы

Построим в «Simulink» непрерывную модель системы управления и произведем настройку непрерывного ПИД-регулятора. Настройку параметров ПИД-регулятора будем производить по кривым управляемого и возмущающего процессов, добиваясь их оптимального вида:

x(t) = 1(t) => y(t) = h(t)

Рис 5.1.2. Схема настройки ПИД регулятора по управляемому процессу

n(t) = 1(t) => y(t) = y_в(t)



Рис 4.1.3. Схема настройки ПИД регулятора по возмущающему процессу

Настройкой параметров ПИД-регулятора добьемся следующих наиболее оптимальных характеристик h(t) и y_в(t):

Рис 5.1.4. Оптимальная характеристика h(t)

Рис 4.1.5. Оптимальная характеристика y_в (t)

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

K_p = 0,1; T_и = 0,1; T_д = 3;

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

Исследуем теперь дискретную модель системы автоматического управления с ПИД-регулятором.

4.2 Настройка ПИД-регулятора дискретной модели системы

Рис 4.2.1. Схема дискретной модели системы

Исследуем дискретную модель системы с ПИД-регулятором, параметры которого попытаемся получить пересчетом параметров непрерывного регулятора из п. 4.1.

q₀ = K_p(1 + ) = 0,1(1+ 30) = 3,1;

q₁ = - K_p(1 - + ) = -0,1(1 - 1 + 60) = - 6;

q₂ = K_p= 0,1·30 = 3;

Построим в «Simulink» модель цифровой системы с рассчитанным регулятором и построим кривые управляемого и возмущающего процессов:

Рис 4.2.2. Схема настройки ПИД регулятора по управляемому процессу

Рис 4.2.3. Схема настройки ПИД регулятора по возмущаемому процессу

При рассчитанных параметрах цифрового ПИД-регулятора получаем следующие характеристики h(k) и y_в(k):

Рис 4.2.4. Характеристика h(k) системы с цифровым регулятором, рассчитанным по параметрам непрерывного



Рис 4.2.5. Характеристика y_в(k) системы с цифровым регулятором, рассчитанным по параметрам непрерывного

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

5. Настройка системы с цифровым регулятором с минимальной обобщенной дисперсией

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

G_р (z)= ;

G_о (z) = = ;

G_фф (z) = = + =

= =

==;

=> G_р(z)= =

==

==

=;

6. Сравнение полученных цифровых систем с разными регуляторами

6.1 Характеристики по возмущающему воздействию

Сравним выходные характеристики систем по возмущающему воздействию:

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

Рис 6.1.2. Выходная характеристика по возмущающему воздействию системы с ПИД-регулятором



Рис 6.1.3. Выходная характеристика по возмущающему

воздействию системы с регулятором с минимальной обобщенной дисперсией

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

6.2 Оценка по минимуму квадратичной интегральной оценки

Оценку по данному критерию произведем, построив в «Simulink» все три модели с использованием блока «Display»:

Рис 6.2.1. Исходная система



Рис 6.2.2. Система с ПИД-регулятором

Рис 6.2.2. Система с регулятором с минимальной обобщенной дисперсией

Из рисунков видно, что лучшей по данному критерию является система с ПИД-регулятором.

7. Рекомендации по выбору регулятора

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

Список литературы

Чостковский Б.К. Методическое пособие. Имитационное моделирование оптимального управления стохастическим объектом.- Самара: СамГТУ.

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