авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 | 2 || 4 |

«Министерство общего и профессионального образования Российской Федерации Тверской государственный технический университет В.Ф. ...»

-- [ Страница 3 ] --

В качестве А-регулятора используется регулятор с дробно рациональной передаточной функцией. Для того, чтобы число уравнений для определения параметров настройки регулятора было не меньше числа параметров регулятора, порядок передаточной функции А-регулятора должен быть равен порядку приведенной непрерывной части. Если, кроме того, мы хотим иметь астатический регулятор, обеспечивающий нулевую статическую ошибку в замкнутой системе, знаменатель передаточной функции регулятора должен содержать сомножитель (z-1) (как и все астатические регуляторы: И-, ПИ-, ПИД) Отметим, что цифровая АСР с А-регулятором имеет конечный, апериодический и оптимальный по длительности переходной процесс только при ступенчатом воздействии.

К недостаткам систем с А-регуляторами можно отнести следующие:

1. Передаточная функция А-регулятора зависит от передаточной функции приведенной непрерывной части и усложняется с ростом порядка последней.

2. С уменьшением периода квантования возрастают максимальное значение регулирующего воздействия и значения параметров настройки регулятора.

В качестве примера рассмотрим цифровую АСР с объектом регулирования в виде инерционного звена первого порядка и ПИ регулятором.

Напомним (см. раздел 6.3), что передаточная функция приведенной непрерывной части, соответствующая данному объекту, равна d WПНЧ ( z ) =, (200) z z где d = K (1 z1 );

z1 = e T Tоб, K, Tоб – параметры инерционного звена первого порядка.

Тогда с учетом передаточной функции ПИ-регулятора (176) пердаточная функция замкнутой системы равна:

d z z Wзс ( z ) = d b1 z + b 1+ z z1 z Характеристическое уравнение замкнутой системы:

( z z1 )( z 1) + d (b1 z + b0 ) = z 2 + z ( z1 1 + db1 ) + ( z1 + db0 ) = Условия конечной длительности переходного процесса (199):

z1 1 + db1 =, z1 + db0 = 0 или с учетом обозначений для b1 и b0 :

d ( K1* + K 0*T ) = 1 + z.

dK1* = z1 Решая последнюю систему уравнений, получаем настройки ПИ регулятора, обеспечивающие конечную длительность переходного процесса:

e T T z1 об K= = *, K (1 z1 ) K (1 e T 1 Tоб ) 1 K 0* = =.

TK (1 z1 ) TK (1 e T Tоб ) При T 0 K1*, K 0*.

Желаемая передаточная функция замкнутой системы (197):

Wжел ( z ) = d ( z 1 z 2 ), т.е. переходной процесс в рассматриваемой системе действительно заканчивается за 2T.

Расчет настроек цифровых регуляторов методом расширенных частотных характеристик Нахождение расширенных частотных характеристик Для нахождения расширенных частотных характеристик следует в выражение для передаточной функции подставить z = e pT = e ( m + j )T = e mT (cos T + j sin T ).

Обозначив для упрощения записи e mT = e;

cos T = c1 ;

sin T = s1, получаем z = e(c1 + s1 ). (201) Приведем в качестве примеров выражения для расширенных частотных характеристик приведенной непрерывной части (200) и ПИ регулятора (176).

РЧХ приведенной непрерывной части (200):

d A* ( ) =, e 2 2ec1 z1 + z es * ( ) = arctg.

ec1 z РЧХ ПИ-регулятора (176):

b12 e 2 + 2b1b0 ec1 + b A ( ) = *, e 2 2ec1 + b1es1 es * ( ) = arctg arctg.

b1ec1 + b0 ec1 Нахождение ЛРЗ Для нахождения уравнения линии равного затухания используются дискретные аналоги условий заданной колебательности переходного процесса в замкнутой системе (36):

AПНЧ (m, ) A* (m, ) = 1 * рег (202) ПНЧ (m, ) + * (m, ) = * рег Поступая так же, как при выводе уравнений ЛРЗ для непрерывных систем, можно получить дискретные аналоги этих уравнений.

Дискретный аналог выражений (41) для П-регулятора:

K1* = (m, 1 ) * A (203) ПНЧ 1 = arg[ ПНЧ (m, 1 )] = * Дискретный аналог выражений (43) для И-регулятора:

e 2 2ec1 + K = * TeAПНЧ (m, ) = 0 * (204) * es 0 = arg ПНЧ (m, ) + T arctg = ec1 1 = Дискретный аналог уравнений (47), (50) ЛРЗ для ПИ-регулятора:

es ПИ (m, ) = ПНЧ (m, ) + arctg * * ec1 1 e 2 2ec1 + 1 K ( ) = sin( ПИ T ) * * (205) s1 AПНЧ (m, ) 1 * e 2 2ec1 + [sin ПИ e sin( ПИ T )] K ( ) = * * * Tes1 AПНЧ (m, ) 0 * Дискретный аналог уравнений (59) для ПИД-регулятора:

e 2 2ec1 + 1 K* c sin ( ПИ T ) + 2 2 1 K ( ) = * * s1 AПНЧ (m, ) T e * (206) e 2 2ec1 + 1 c1 K 2* K ( ) = sin ПИ e sin( ПИ T ) + 2 1 2 + * * * Tes1 AПНЧ (m, ) T e e * Особые прямые Специфической особенностью ЛРЗ для дискретных систем является то, что при = T ЛРЗ вырождается в особую прямую (рис. 74) K 0* b T 0 K1* Рис. 74.

Запишем уравнение особой прямой для ЛРЗ для системы с ПИД регулятором.

При = T e = e m, c1 = 1, s1 = 0, e 2 2ec1 + 1 = 1 + e m, * ПИ m, = ПНЧ m, = и уравнения (206) вырождаются в особую * T T прямую с уравнением:

* * K 0 = aK1 + b, T T a = 1 T (207) 1 K 2* 1 + e 2m 1 + e m b= + T m * T e 2m e AПНЧ m, T (На рис. 74 tg = a ).

Следует заметить, что если 1 T, т.е. период квантования достаточно мал и выбран правильно, частота T оказывается за пределами частотного диапазона, в котором строится ЛРЗ, поэтому наличие особой прямой не сказывается на форме ЛРЗ. Если же период квантования слишком велик и выбран неправильно, т.е. 1 T, ЛРЗ вырождается в особую прямую.

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

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

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

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

7.1. Основные характеристики случайных процессов [8, 9, 7] Случайной функцией времени (случайным, стохастическим процессом) x(t) называется функция, в любой момент времени представляющая случайную величину. Дискретным аналогом случайного процесса является временная последовательность (временной ряд) xi = x(iT ), i = 0,1,2, K, где Т – по-прежнему период отсчетов сигнала (интервал дискретности, период квантования).

Если случайная величина x(t0) в произвольный момент времени t (или произвольный член временного ряда xi ) имеет нормальное распределение, то статистические свойства случайного процесса (временного ряда) полностью определяют четыре характеристики:

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

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

Математическое ожидание Математическое ожидание случайного процесса характеризует его истинное среднее значение:

M x = M {x(t )} (M – символ операции нахождения математического ожидания) Математическое ожидание временного ряда:

M x = M {xi } Оптимальной оценкой математического ожидания стационарного эргодического временного ряда является среднее арифметическое отсчетов:

1 N € x x = Mx =, (208) i N i = N – число членов временного ряда.

Рекуррентная оценка среднего арифметического:

xn = xn1 + ( xn xn1 ), (209) N xn, xn1 – значения среднего арифметического в текущий и предыдущий моменты времени.

Дисперсия Дисперсия случайного процесса (временного ряда) – это математическое ожидание квадрата отклонения случайного процесса (временного ряда) от его математического ожидания. Другими словами, дисперсия – это средний квадрат отклонения случайного процесса (временного ряда) от его математического ожидания:

{ } Dx = x2 = M [x(t ) M x ], (210) или для временного ряда { } Dx = M [xi M x ] Оптимальной оценкой дисперсии временного ряда является среднее арифметическое квадратов отклонений измеренных значений от оценки математического ожидания:

1N 1N € ( xi x ) = N 1 ( xi x 2 ) s = Dx = 2 (211) N 1 i x i Рекуррентная формула для оценки дисперсии:

N 2 2 s n1 + ( xn xn1 ) s n2 = N 1 N Итак, математическое ожидание характеризует истинное среднее значение случайного процесса (временного ряда), а дисперсия – его случайный разброс относительно математического ожидания.

Поскольку дисперсия имеет размерность квадрата случайного процесса, вводят дополнительную характеристику – среднеквадратическое отклонение (СКО) случайного процесса, определяемое как положительное значение корня квадратного из дисперсии:

x = Dx.

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

Rx ( ) = M {x(t + ) x(t )} (212) Данное определение предложено Н. Винером (инженерное определение корреляционной функции).

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

C x ( ) = M {[x(t + ) M x ][x(t ) M x ]} (213) Ковариационная функция связана с корреляционной через квадрат математического ожидания случайного процесса:

C x ( ) = Rx ( ) M x2 (214) Для центрированных процессов, математическое ожидание которых равно нулю:

M [x(t ) M x ] = 0, ковариационная функция равна корреляционной.

Сопоставляя (210), (213), (214), получаем C x (0) = Dx = Rx (0) M x2, (215) или Rx (0) = Dx + M x2.

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

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

Rx ( ) = Rx ( ), т.к. знак (направление) сдвига не имеет значения, поэтому их можно строить только при 0. Кроме того, корреляционная функция стационарного процесса является функцией, убывающей от Dx + M x2 до M x2 и имеющей при этом монотонный или колебательный характер.

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

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

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

Rx ( ) 4 Рис. 75.

Линиям 2 и 3 на рис. 75 соответствуют корреляционные функции низкочастотного и высокочастотного процессов.

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

Rx (kT ) = M [xi + k xi ], где k – число периодов квантования во временном сдвиге (лаге) :

= kT.

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

1 N k € x Rx ( k ) = xi, (216) i+k N k i = где N – число членов временного ряда.

С ростом лага k точность вычисления оценки (216) падает, т.к.

уменьшается число членов ряда, по которому считается оценка. При k= € оценка Rx (0) считается как среднее арифметическое N отсчетов, а при k=N-2 среднее арифметическое считается всего по двум точкам, т.е. весьма неточно. Поэтому для того, чтобы оценка корреляционной функции была достаточно точной, должно выполнятся условие kmax N, например, kmax=(0.01–0.05)N.

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

Спектральная плотность мощности S x ( ) есть предел отношения мощности N случайного процесса на выходе узкополосного фильтра с полосой к величине этой полосы при 0 :

N S x ( ) = lim Таким образом, физический смысл спектра мощности – это скорость изменения (плотность распределения) мощности сигнала по частоте.

Спектр мощности есть неотрицательная функция частоты: S x ( ) 0.

Чем более низкочастотным является случайный сигнал, тем быстрее затухает его спектральная плотность. На рис. 76 кривая спектральной плотности 2 соответствует более низкочастотному сигналу, чем кривая 3.

S() 1 Рис. 76.

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

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

Коррелограммный способ (Преобразование Фурье корреляционной функции) В англоязычной литературе корреляционная функция называется коррелограммной, что и объясняет название способа.

Спектральная плотность мощности и корреляционная функция связаны парой прямого и обратного преобразования Фурье:

S x ( ) = Rx ( )e j d S ( )e d Rx = j x Экспериментальная оценка спектра мощности рассчитывается на интервале конечной длительности t p, поэтому бесконечные пределы интегрирования в формулах следует заменить конечными:

tp € € S x ( ) = Rx ( )e j d (217) t p С учетом формулы Эйлера:

e j = cos j sin, получаем для (217):

tp 2 tp € € € S x ( ) = Rx ( ) cos d j Rx ( ) sin d.

(218) t p 2 t p Второй интеграл в правой части (218) равен нулю, как интеграл от нечетной функции в симметричных пределах, следовательно, tp € € S x ( ) = Rx ( ) cos d, t p или для вещественных данных tp € € S x ( ) = 2 Rx ( ) cos d.

При вычислении оценки спектральной плотности мощности на ЦВМ дискретным аналогом формулы (217) является выражение L € € w(i) R (i)e S x ( ) = jiT, x i = L где L – максимальный лаг при вычислении оценки корреляционной функции;

w(i) – корреляционное окно – специальная функция, с помощью которой осуществляется «взвешивание» ординат корреляционной функции с целью улучшения оценки спектральной плотности мощности. В частном случае при w(i) = 1 имеем прямоугольное корреляционное окно.

Считая, что = k, k = 0, N (k – номер ординаты оценки спектра мощности, N – количество отсчетов временного ряда), окончательно получаем следующую формулу для вычисления оценки спектральной плотности мощности:

L j ik € € S x (k ) = w(i ) Rx (i )e N, k = 0, N 1 (219) i = L Определение спектральной плотности мощности преобразованием Фурье параметрической модели временного ряда Многие встречающиеся на практике временные ряды удовлетворительно описываются следующим конечно-разностным уравнением, называемым моделью авторегрессии – скользящего среднего (АРСС):

xn + a1 xn1 + K + a p xn p = b0u n + b1u n1 + K + bq u nq, p q или xn = ai xni + bi u k i (220) i =1 i = Первая сумма в правой части (220) называется авторегрессионной частью модели АРСС и характеризует зависимость временного ряда в текущий момент xn от прошлых значений временного ряда xn-1, взвешенных с коэффициентами ai.

p – порядок авторегрессии (АР).

Вторая сумма в правой части (220) соответствует модели скользящего среднего.

un – белый шум, возбуждающий модель АРСС, с нулевым математическим ожиданием и дисперсией u2.

q – порядок модели скользящего среднего (СС).

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

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

Взяв Z-преобразование модели АРСС (220), можем найти дискретную передаточную функцию формирующего фильтра:

x( z ) B( z ) Wф ( z ) = =, u ( z ) A( z ) где q B( z ) = 1 + bi z i i =, p A( z ) = 1 + ai z i i = полиномы от z-i порядков q и p соответственно.

Частными случаями модели АРСС являются модели авторегрессии (q = 0) и скользящего среднего (p = 0).

Учитывая, что, как известно из теории управления, спектры мощности выходного и входного сигналов динамического звена связаны через квадрат его АЧХ, а также, что спектр мощности дискретного белого шума на входе формирующего фильтра равен S u ( ) = u2T, получаем для спектральной плотности мощности временного ряда, описываемого моделью АРСС (220), следующее выражение:

B ( j ) B ( j ) B ( j ) S x ( ) = T = u2T, (221) A( j ) A( j ) A( j ) u где A( j ), B ( j ) - изображения Фурье полиномов A( z ), B( z ).

7.2. Определение дисперсии выходной величины в цифровой АСР [10, 11] Считаем, что возмущающее воздействие в цифровой системе представляет случайную временную последовательность, описываемую моделью АРСС (220). Спектральная плотность мощности этого воздействия определяется соотношением (221).

Для удобства дальнейших выводов перейдем из частотной плоскости в плоскость z. Аналогом соотношения (221) в плоскости z является выражение:

B ( z ) B ( z 1 ) S x ( z) = T (222) A( z ) A( z 1 ) u (Напомним, что такой переход осуществляется с помощью подстановки z = e jT ).

Например, для модели авторегрессии первого порядка АР1 (p = 1, q = 0):

xn = a1 xn1 + u n, имеем:

B ( j ) = 1;

A( j ) = 1 + a1e jT, A( z ) = 1 + a1 z 1, A( z 1 ) = 1 + a1 z.

Тогда согласно (221) и (222) u2T S x ( ) =, (223) 1 + 2a1 cos T + a u2T S x ( z) = 1 + 2a1 ( z + z 1 ) + a Как отмечалось, спектральные плотности выходной и входной случайных последовательностей динамической системы связаны через квадрат её АЧХ. Переходя в плоскость z выражение для спектра мощности выхода системы, можно записать следующим образом:

S y ( z ) = Wзс ( z ) S x ( z ) = Wзс ( z )Wзс ( z 1 ) S x ( z ), (224) где S y ( z ), S x ( z ) - спектральные плотности временных рядов на выходе и входе системы, Wзс (z ) - дискретная передаточная функция замкнутой системы.

В разделе 7.1. также отмечалось, что корреляционная функция и спектр мощности связаны парой преобразований Фурье (или при переходе в выражении (222) в плоскость z – парой Z-преобразования). Поэтому корреляционная функция выходной величины цифровой системы может быть найдена с помощью обратного Z-преобразования спектральной плотности выходной величины, которое определяется следующим образом:

S ( z ) z k 1 dz R y (k ) = 2j y z = или с учетом (224) W Ry ( k ) = ( z ) S x ( z ) z k 1dz, (225) зс z = - круговой интеграл по контуру z = 1, т.е. интеграл, в котором где интегрирование осуществляется по замкнутому контуру, т.к. при изменении частоты от 0 до 2 T конец вектора z = e jT описывает в плоскости z окружность с единичным радиусом.

Согласно (215) дисперсия выходной величины равна y2 = Ry (0) M y или с учетом (225) окончательно получаем:

y2 = W ( z ) S x ( z ) z 1dz M y2 (226) 2j зс z = (Ниже для простоты считаем M y = 0 ).

Для нахождения интеграла (226) можно использовать два способа:

численное интегрирование и использование итерационной процедуры.

Численное вычисление интеграла (226) в частотной области Для перехода в частотную область опять воспользуемся подстановкой z = e jT.

Тогда учитывая, что dz = e jT jTd и z 1dz = jTd, а также то, что в силу периодичности подинтегральной функции интеграл в пределах периода 0 2 T равен удвоенному интегралу в пределах полупериода 0 T, получаем:

T T y2 = Wзс ( j ) S x ( )d.

* (227) В частности, при численном интегрировании методом прямоугольников, обозначая = i ;

T = N, откуда = NT ( – шаг квантования по частоте, N – число интервалов дискретности), сводим интеграл (227) к сумме:

1 N W y2 = ( ji ) S x (i ).

* (228) зс N i = В качестве примера запишем выражение (228) для системы второго порядка, на входе которой действует случайная последовательность, описываемая моделью АР1.

При = i, = NT, T = i N, поэтому модель (223) принимает форму u2 T S x (i ) = (229) 1 + 2a1 cos(i N ) + a Пусть c1 z + c Wзс ( z ) =, d 2 z + d1 z + d C ( z ) = c1 z + c0, D ( z ) = d 2 z 2 + d1 z + d 0.

Тогда D ( z ) = D ( z ) D ( z 1 ) = (d 02 + d12 + d 22 ) + (d 0 d1 + d1d 2 )( z + z 1 ) + d 0 d 2 ( z 2 + z 2 ) Переходя в последнем выражении с помощью подстановки z k + z k = 2 cos(kT ) = 2 cos(ki N ) в частотную область, получаем D ( ji ) = (d 02 + d12 + d 22 ) + 2(d 0 d1 + d1d 2 ) cos(i N ) + 2d 0 d 2 cos(2i N ).

(230) Действуя аналогично, находим C ( ji ) = (c02 + c12 ) + 2c0 c1 cos(i N ) (231) и D( ji ) W ( ji ) = * (232) C ( ji ) зс В общем случае для полинома p-го порядка p D( z ) = d k z k k = квадрат его модуля можно определить по выражению:

p 1 p p D( ji ) = d k2 + 2( d k d k +1 ) cos(i N ) + 2( d k d k + 2 ) cos(2i N ) + K k =0 k =0 k = p k p p + 2d 0 d p cos( pi N ) = d k2 + 2 cos(ki N ) d s d s + k k =0 k =1 S = Подставляя, (230) и (231) в (232), а (229) и (232) в (228), получаем искомое выражение для дисперсии выходной величины y2.

Вычисление интеграла (226) с использованием итеративной процедуры Определение дисперсии y2 по выражению (226) сводится к решению интеграла вида:

B( z ) B ( z 1 ) I= z dz, (233) 2j z =1 A( z ) A( z ) где полиномы A(z ) и B(z ) имеют вид A( z ) = a0 z n + a1 z n1 + K + an B( z ) = b0 z n + b1 z n1 + K + b Обозначим A* ( z ) = z n A( z 1 ), (234) Ak ( z ) = a0k z k + a1k z k 1 + K + akk (235) Bk ( z ) = b0k z k + b1k z k 1 + K + bkk An ( z ) = A( z );

Bn ( z ) = B( z ) Полиномы Ak ( z ), Bk ( z ) определяются с помощью следующих рекуррентных соотношений:

Ak 1 ( z ) = z 1 [Ak ( z ) k Ak* ( z )] (236) Bk 1 ( z ) = z 1 [Bk ( z ) k Ak* ( z )];

k = 1, n, где bkk akk k = k ;

k = k.

a0 a Из выражений (234)(236) следует, что коэффициенты многочленов Ak ( z ), Bk ( z ) задаются следующими рекуррентными формулами:

aik 1 = aik k akki (237) bik 1 = bik k akki, i = 0, k с начальными условиями ain = ai ;

bin = bi.

Уравнения (237) имеют смысл, если a0k 0. Коэффициент a0n всегда можно выбрать ненулевым.

Интеграл (233) можно вычислить по рекуррентным формулам, используя последовательность интегралов:

Bk ( z ) Bk ( z 1 ) Ik = z dz, k = 0, n, (238) 2j z =1 Ak ( z ) Ak ( z ) In = I Рекуррентная формула для определения I k имеет вид I k = K + I k 1 (1 K ), 2 I 0 = 02.

Используя формулы 237, можно получить следующее выражение для интеграла (238):

(bii ) 1 k ai Ik = (239) a0k i =0 Для удобства вычислений по формуле (239) можно использовать таблицу 8.

Таблица Коэффициенты ain Коэффициенты bin a0 a1.... an1 an b0 b1.... bn1 bn * an an1 a1 a0 an an1 a1 a a0n1 a1n1 b0n1 b1n1 bnn ann * a0n1 a0n ann11 ann1 ann11 ann 2 M M a0 a b0 b * a11 a a11 a a00 b * В табл. 8 кружками помечены коэффициенты, входящие в выражение (239).

k Каждая четная строка коэффициентов ai получается путем записи коэффициентов предыдущей строки (отмечены знаком ) в обратном k k порядке. Четные строки коэффициентов ai и bi совпадают. Элементы k k нечетных строк коэффициентов ai и bi вычисляются с помощью соотношений (237) 7.3. Синтез регулятора с минимальной дисперсией [5] Целью расчета регуляторов с минимальной дисперсией (МД – регуляторов) является минимизация дисперсии выходной величины цифровой системы или равного ей с точностью до M y математического ожидания квадрата выходной величины:

{} min y = min M y k 2 Поскольку решение этой задачи имеет смысл только при ограничениях на управления (иначе, увеличивая управляющее воздействие, дисперсию выходной величины можно сделать сколь угодно малой) в критерий управления I следует ввести ещё одно слагаемое – средний квадрат управления с весовым коэффициентом r :

{ } I = M y k +1 + rxk min 2 (240) В критерии (240) используется последующее значение выхода y k + (а не текущее y k ), т.к. текущее значение управления x k определяет в силу запаздывания в объекте управления последующее значение выхода y k +1.

Рассмотрим в качестве примера структурную схему цифровой АСР со случайным воздействием F, приведенным к выходу приведенной непрерывной части (рис. 77) F u(z) Wф (z ) y зад e(z ) y (z ) x(z ) W рег (z ) W ПНЧ (z ) – Рис. 77.

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

b0 z 1 B( z ) WПНЧ ( z ) = =.

a1 + a0 z 1 A( z ) B( z ) = b0 z 1, A( z ) = a1 + a 0 z 1, а случайное воздействие описывается моделью АР1:

Fk = cFk 1 + u k.

Передаточная функция формирующего фильтра для этой модели:

F ( z) 1 Wф ( z ) = = =, u ( z ) 1 + cz C ( z) C ( z ) = 1 + cz 1.

Передаточную функцию МД–регулятора будем также искать в виде дробно–рациональной функции от z :

x( z ) W рег ( z ) =, e( z ) или, считая y зад = 0 и ek = y зад y k = y k, x( z ) W рег ( z ) = (241) y( z) Для нахождения выражения для выходной величины y k +1, входящей в критерий (240), запишем её Z – изображение:

zy ( z ) = WПНЧ ( z ) zx( z ) + Wф ( z ) zu ( z ) = B( z ) = zx( z ) + zu ( z ) (242) A( z ) C ( z) Подставляя в (242) выражения для полиномов A( z ), B( z ), C ( z ), приводя результаты к общему знаменателю и совершая обратное Z – преобразование, получаем:

a1 y k +1 + (a 0 + a1c ) y k + a 0 cy k 1 = b0 x k + b0 cx k 1 + a1u k +1 + a 0 u k (243) Определяя из последнего выражения y k +1 и подставляя его в критерий (240), имеем:

a +ac a a0 c b0 b0 c I = M yk y k 1 + xk + xk 1 + u k +1 + u k + 0 a1 a1 a1 a1 a + rxk min (244) В текущий момент времени k известны текущее u k и прошлые значения возбуждающего шума u k i. Неизвестно будущее значение u k +1.

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

M [u k i ] = u k i, i 0, т.к. математическое ожидание неслучайной величины равно самой этой величине. В то же время математическое ожидание последующего значения шума M [u k +1 ] равно нулю, поскольку возбуждающий шум по предположению имеет нулевое математическое ожидание.

a – константа, а x – случайная величина с нулевым Если математическим ожиданием, имеет место очевидное соотношение:

M [(a + x ) ] = a 2 + M [ x 2 ] Поэтому после нахождения математического ожидания выражения в фигурных скобках (244) получаем:

a + a1c ac b bc a I = 0 y k 0 yk 1 + 0 xk + 0 xk 1 + 0 u k + a1 a1 a1 a1 a + M [u k +1 ] + rxk min 2 Оптимальное значение управления xk определяем из необходимого условия экстремума критерия оптимальности:

a +a c b I ac b bc a = 2 0 1 y k 0 y k 1 + 0 xk + 0 xk 1 + 0 u k 0 + 2rxk = 0 (245) xk a1 a1 a1 a1 a1 a Для нахождения передаточной функции МД–регулятора запишем Z – преобразование соотношения (245). Выражение в квадратных скобках (245) есть, как следует из (243), y k +1 за вычетом u k +1, поэтому его Z – преобразование равно:

zy ( z ) zu ( z ), следовательно, Z – преобразование (245) принимает вид [zy( z ) zu ( z )] b0 + rx( z ) = 0 (246) a Найдем из (242) zu (z ) и подставим в (246).

B ( z )C ( z ) zu ( z ) = C ( z ) zy ( z ) zx( z ) (247) A( z ) После подстановки (247) в (246), имеем b B( z ) zy ( z ) C ( z ) zy ( z ) + C ( z ) A( z ) zx( z ) a + rx( z ) = Передаточную функцию МД–регулятора получим разрешив последнее равенство относительно отношения x( z ) y ( z ), как того требует определение (241):

zA( z )[C ( z ) 1] b0 a W рег ( z ) =, rA( z ) zC ( z ) B( z ) b0 a или с учетом выражений для полиномов A( z ), B( z ), C ( z ) :

(a ) + a0 z 1 cb0 a W рег ( z ) = ( )( ) (248) ra1 b02 a1 + ra0 cb02 a1 z 8. Синтез многомерных дискретных регуляторов в пространстве состояния [12, 13, 14] 8.1. Формулировка задачи оптимального управления Задача синтеза многомерного регулятора в пространстве состояния формализуется обычно в виде задачи о максимальной точности с интегральным квадратичным функционалом I :

t ( ) r rrr min J = x T Q x + u T R u dt, (249) u t где t 0, t1 – соответственно начальное и конечное время, (t1 t 0 ) – интервал управления (при t1 имеем бесконечный интервал управления), r x – вектор переменных состояния системы размерностью n ( n – вектор);

r u – r – вектор переменных управления;

Q – ( n n ) – квадратная положительно определенная диагональная матрица весовых коэффициентов состояния;

R – ( r r ) – квадратная положительно полуопределенная диагональная матрица весовых коэффициентов управления.

r Вообще говоря, в критерии задачи (249) вместо вектора состояния x следовало бы писать отклонение вектора состояния от его заданного r r значения x зад. Однако для упрощения записи полагаем x зад = 0.

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

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

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

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

( ) N r r rr min I = xiT 1Q xi +1 + uiT R ui, (250) + i = где N – число периодов квантования по времени в интервале управления (Как в задаче (240), последующее состояние и предыдущее управление в силу запаздывания в объекте управления отнесены к одному шагу).

8.2. Уравнения состояния и измерения При решении задачи оптимального управления (249) уравнения динамики объекта (системы) управления записывают в виде системы дифференциальных уравнений первого порядка в матричной форме записи и называют матричным уравнением состояния.

Для линейных стационарных систем уравнение состояния имеет вид:

r r r & x = Ax + Bu, (251) r & где x – n – вектор первых производных переменных состояния, A и B – соответственно ( n n ) – и ( n r ) – матрицы состояния и управления.

В реальных условиях измерению доступны только отдельные элементы вектора состояния или их линейные комбинации, поэтому уравнение состояния дополняется матричным уравнением измерения r (выхода), связывающим l – вектор выходных переменных y с вектором состояния:

r r y =Cx, (252) где C – ( l n ) – матрица выхода.

В качестве примера рассмотрим получение матричного уравнения состояния для одномерного по управлению объекта, описываемого дифференциальным уравнением n – го порядка (2), которое перепишем в виде d n 1 y d m 1u dny d mu a n n + a n 1 n 1 + L + a 0 y = bm m + bm 1 m 1 + L + b0 u (253) dt dt dt dt Пусть для определенности m = n 1.

Запишем преобразование Лапласа уравнения (253) (a ) ( ) p n + a n 1 p n 1 + L + a 0 y ( p) = bn 1 p n 1 + bn 2 p n 2 + L + b0 u ( p), (254) n или A( p) y ( p) = B( p)u ( p ), (255) p порядков n и ( n 1 ) A( p) и B( p) – полиномы от где соответственно.

Из (255) следует, что y ( p) = B( p) u ( p) A( p), (256) Обозначим u ( p) A( p) = x( p), (257) тогда выражения (256), (257) можно записать в виде:

y ( p) = B( p ) x( p) (258) A( p) x( p) = u ( p) (259) Запишем оригиналы выражений (258), (259) n 1 n y = bn1 x + bn2 x + L + b1 x + b0 x & (260) n n a n x + a n 1 x + L + a1 x + a 0 x = u & (261) Введем переменные состояния x1 = x x2 = x& x3 = && x LL (262) n xn1 = x n xn = x Тогда уравнение (261) с учетом обозначений (262) можно записать в виде следующей системы уравнений первого порядка:

x1 = x2 & x2 = x & LL (263) xn1 = xn & (u a0 x1 a1 x2 L an1 xn ) xn = & an По существу первые ( n 1 ) уравнений системы (263) являются обозначениями, а n – ое уравнение получено из уравнения (261).

Вводя обозначения 0 1 0 0 L 0 0 1 0 L =L L L L L =M Ann Bn ;

, 0 0 0 L M a a a a 0 1 2 n L an a1 an a0 an получаем матричное уравнение состояния в форме (251).

A имеет в данном случае так называемую форму Матрица Фробениуса.

Уравнение (260) с учетом обозначений (262) можно записать в виде:

y = b0 x1 + b1 x 2 + L + bn 1 x n. (264) Обозначая C1n = b0 b1 L L bn2, получаем уравнение выхода в форме (252). В данном случае y – скалярная переменная, являющаяся линейной комбинацией переменных состояния.

Матричная форма записи уравнений состояния – выхода является удобным инструментом при анализе многосвязных систем. Например, для двумерного объекта, структурная схема которого изображена на рис. ( n = 5, r = 2, l = 2 ), u1 x1 x2 y K1 K T1 p + 1 p K x K T3 p + u2 x4 x6 y K4 K T4 p + 1 T5 p + Рис. 78.

описываемого следующими системами уравнений:

x1 = (1 T1 )x1 + (K 1 T1 )u & x 2 = K 2 x & x3 = (1 T3 )x3 + (K 3 T3 )x &, & 4 = (1 T4 )x 4 + (K 4 T4 )u x x5 = (1 T5 )x5 + (K 5 T5 )x & y1 = x 2 K 6 x y 2 = x3 + x5, матрицы состояния, управления и выхода выглядят следующим образом:

1 T4 0 0 0 0 K 1 T1 0 0 K2 0 0 0 0 0 0 1 A = K 3 T3, B=, CT = 0 1 T3 0 0 0 0 1 T4 K 0 0 0 0 0 K 4 T4 1 T 0 0 0 K 5 T5 0 0 0 Для перехода к дискретному уравнению состояния рассмотрим решение уравнения состояния (251) при постоянных матрицах A, B, r r начальном состоянии x0 и ступенчатом управлении u 0 :

[ ] r r r x (t ) = e A(t t0 ) x0 + A1 e A(t t0 ) E Bu0 (265) e A(t t0 ) Матричная экспонента определяется разложением в степенной ряд e A(t t0 ) E + A(t t 0 ) + A2 (t t 0 ) 2 2! + L + Ak (t t 0 ) k k ! + L (266) Обозначим Ft t0 = e A( t t0 ) (267) [ ] Gt t0 = A1 e A(t t0 ) E B = A1 (F E )B (268) Матрицы F и G называются переходными матрицами состояния и управления соответственно.

Если матрица A плохо обусловлена, нахождение матрицы B вызывает определенные трудности, поэтому желательно определить матрицу G так, чтобы избежать операции обращения матрицы A.

Введем в рассмотрение матрицу :

( ) t t0 = A1 Ft t0 E С учетом (266), (267) выражение для матрицы принимает вид:

t t0 = E (t t 0 ) + A(t t 0 ) 2 2 ! + L + Ak 1 (t t 0 ) k k ! + L (269) Тогда матрицы G и F могут быть выражены через матрицу следующим образом:

Gt t0 = t t0 B (270) Ft t0 = At t0 + E (271) Как видим, для вычисления матрицы G по формуле (270) не требуется обращать матрицу A.

С учетом обозначений (267), (268) выражение (265) для переходной функции приобретает вид:

r r r x (t ) = Ft t0 x0 + Gt t0 u Положив в последнем соотношении r r r r t 0 = kT ;

x0 = x[kT ] ;

u0 = u[kT ] ;

t = (k + 1)T ;

t t 0 = T, получаем дискретное уравнение состояния:

r r r x[(k + 1)T ] = F (T ) x[kT ] + G (T )u[kT ] или в упрощенной форме записи:

r r r x k +1 = Fx k + Gu k (272) Матрицы F (T ) и G (T ) можно найти из выражений (269) (271), заменив в них (t t 0 ) на T :

(T ) = ET + AT 2 2 ! + A2T 3 3! + L + A k 1T k k ! + L (273) F (T ) = A (T ) + E (274) G (T ) = (T ) B (275) Дискретным аналогом уравнения выхода (252) является уравнение r r y k = Cx k (276) 8.3. Синтез дискретного П–регулятора состояния Решаем задачу (250) при условии, что объект управления описывается уравнением (272).

Предположим, что в процессе решения мы находимся на k –том от начала или j –том от конца шаге управления ( k = N j ). Обозначим минимальное значение целевого функционала на оставшихся шагах управления (функцию оптимального поведения) через I k :

(x ) N rT r rT r I = min Q xk +1 + u k R u k (277) k + k uk Lu N k =N j Уравнение дискретного П – регулятора будем искать в виде:

r r u k = K p, k xk, (278) K p, k (r n) матрица где коэффициентов передачи оптимального регулятора на k –том шаге.

r r r r Последовательно выражая x k +1 через x k и u k через x k с помощью уравнений объекта (272) и регулятора (278) для k = N 1, N j, функцию оптимального поведения (277) можно представить в виде квадратичной * формы некоторой матрицы Pk :

rT r I k = xk Pk* xk * (279) или для момента времени ( k + 1 ):

rT r I k +1 = xk +1 Pk*+1 xk + * (280) Pk* – квадратная симметричная ( n n ) – матрица.

В соответствии с принципом оптимальности Беллмана, каково бы ни r было состояние системы перед очередным k –тым шагом, управление u k на этом шаге выбирается так, чтобы минимизировать сумму функции в круглых скобках (250) на k –том шаге и функции оптимального поведения на оставшихся ( k 1 )–м шагах, т.е.

( ) rT r rT r I k = min xk +1Q xk +1 + u k R u k + I k +1, * * u k или с учетом (280):

[( ] ) rT r rT r I k = min xk +1 Q + Pk*+1 xk +1 + u k R u k * (281) uk r Поскольку состояние x k +1 в момент k ещё неизвестно, его можно определить из уравнения объекта (272). При этом выражение (281), в котором для кратности обозначено Pk +1 = Q + Pk*+1, (282) принимает вид:

[ rT r rT r I k = min xk F T Pk +1 F xk + 2 xk F T Pk +1G u k + * uk )] ( rT r + u k G T Pk +1G + R u k (283) Оптимальное управление на k –том шаге находим из условия * стационарности I k :

[ )] I k ( * r r r = 2 G Pk +1 F xk + G Pk +1G + R u k = 0, T T u k откуда ( ) r r u k = G T Pk +1G + R G T Pk +1 F xk (284) Сопоставляя (284) с (278), убеждаемся, что ( ) K p, k = G T Pk +1G + R G T Pk +1 F. (285) Для определения матрицы Pk +1 подставим (284) в (283). После преобразований получаем [ ] ( ) rT r I k = xk F T Pk +1 E G G T Pk +1G + R G T Pk +1 Fxk * (286) Сопоставляя (286) с (279), а также учитывая (282), находим:

[ ]F ( ) Pk = Q + F T Pk +1 E G G T Pk +1G + R G T Pk +1 (287) Уравнение (287) есть конечно–разностный аналог дифференциального уравнения Риккати. При бесконечном интервале управления ( N ) оно вырождается в алгебраическое уравнение:

[ ]F ( ) P = Q + F T P E G G T PG + R G T P (288) а выражение (285) для матрицы коэффициентов передачи регулятора принимает вид:

( ) K p = G T PG + R G T PF (289) Выражения (289), (288) позволяют рассчитать матрицу коэффициентов оптимального по точности П – регулятора состояния.

Уравнение (288) является нелинейным и в общем случае его решение в явном виде получить не удается. Однако его можно рассматривать как рекуррентное соотношение:

Pi = f (Pi 1 ), P0 = E и использовать для его решения итерационные методы. В качестве критерия останова можно использовать, например, близость евклидовых норм матриц на соседних шагах.

Структурная схема многомерной дискретной системы с объектом управления (272) и оптимальным по точности П – регулятором состояния (289) приведена на рис. 79.

r r r xk +1 xk Fk Ee pT G F r T T xзад – (G PG + R ) G T PF T Фиксатор Рис. 79.

8.4. Синтез дискретного ПИ–регулятора состояния – выхода Поскольку система с П–регулятором, как отмечалось в разделе 2.3, характеризуется статической ошибкой, рассмотрим теперь систему с ПИ– регулятором, описываемым уравнением:

k r r r u k = K 0 ei + K1 xk (290) i = rr r ei = y зад yi – l -вектор рассогласования выхода относительно его где r заданного значения (для простоты считаем y зад = 0 );

K1 – соответственно ( r l ) – и ( r n ) – матрицы K0 и коэффициентов передачи И–составляющей и П–оставляющей закона регулирования.

Регулятор, описываемый уравнением (290), является пропорциональным по состоянию и интегральным по выходу и называется поэтому ПИ–регулятором состояния – выхода.

Для того, чтобы избавиться в (290) от суммы, запишем уравнение ПИ–регулятора в приращениях:

r r r rr u k = u k 1 + K 0 ek + K1 ( xk xk 1 ) (291) Будем решать задачу:

( ) N r r rT r min I = ekT+1Qek +1 + u k Ru k, (292) u k k = r r r где u k = u k u k 1 – приращение вектора управления. (Считаем, что объект управления по–прежнему описывается уравнением (272)).

Поскольку, как видно из (290), управление зависит теперь не только от состояния, как в регуляторе (278), но и от суммы рассогласований rp выхода, введем расширенный вектор состояния x k :

k r e i i = r xkp = LL r xk Для того, чтобы избавиться от суммы в определении расширенного вектора состояния, перейдем к его приращению:

r ek r x kp = LL, r x k r r r где x k = x k x k 1.

rp Определим теперь конечно – разностное уравнение для вектора x k.

r Учитывая, что y зад = 0, приращение вектора рассогласования равно:

r r ek +1 = y k +1, или с учетом уравнения выхода (276) r r r ek +1 = ek Cx k +1 (293) Подставляя в (293) конечно – разностное уравнение объекта (272), записанное в приращениях r r r x k +1 = Fx k + Gu k, (294) получаем r r r r ek +1 = ek CFx k CGu k (295) Вводя обозначения E M CF CG F = L L LL ;

G p = LL, p F 0M G можем объединить уравнения (295), (294) в конечно – разностное уравнение эквивалентного объекта (расширенное уравнение состояния):

r r r x kp+1 = F p x kp + G p u k (296) Теперь задачу (292) можно свести к задаче (250), записав её в виде:

((xr ) )Q N r rT r min I = pT xkp+1 + u k Ru k p (297) k + uk k = Поскольку в исходной задаче (292) весовые коэффициенты относились только к вектору рассогласования выхода, расширенная матрица весовых коэффициентов состояния в задаче (297) определяется следующим образом QM Qp = L L L, 0M Таким образом, задача синтеза ПИ–регулятора состояния – выхода (291) сведена к задаче (297) эквивалентной задаче синтеза П–регулятора с расширенным вектором состояния и эквивалентным объектом (296). Решая эту задачу с помощью соотношений, аналогичных (289), (288), получаем r (l + n) – матрицу коэффициентов ПИ–регулятора K ПИ, которую можно расчленить на подматрицы коэффициентов П– и И– составляющих регулятора K 1 и K 0 :

K ПИ = K 0 M K 8.5. Синтез дискретного наблюдателя состояния Поскольку реализация регуляторов (278) или (291) требует знания вектора состояния, возникает задача восстановления состояния по выходным переменным. Эта задача решается с помощью устройств, называемых наблюдателями состояния. Одним из возможных вариантов реализации наблюдателя является наблюдатель Люенбергера, описываемый в дискретном случае уравнением:

[ ] r r r r r € € € xk +1 = Fxk + Gu k + K H y k Cxk, (298) или r r r r xk +1 = (F K H C )xk + Gu k + K H y k, € € (299) r € где x k – оценка вектора состояния, K H (n l) матрица коэффициентов передачи наблюдателя, в определении которой и заключается задача синтеза наблюдателя.

Уравнению (298) соответствует структурная схема на рис. 80.

r uk G r r r € € xk xk + yk Ee pT KH r – € yk F C Рис. 80.

Как видно из рис. 80, наблюдатель состояния представляет r r замкнутую систему с двумя входами y k и u k. Если эта система устойчива и обладает малой статической ошибкой, то по окончании переходного r r r r € € процесса y k стремится к y k и, следовательно, x k стремится к x k.

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

Рассмотрим ещё раз систему с объектом (272) и стационарным регулятором (278). Подставляя (278) в (272), получаем уравнение замкнутой системы:

r r r xk +1 = Fxk + GK p xk, или xk +1 = (F + GK p ) xk, r r которому соответствует следующее характеристическое уравнение:

[ ] det zE ( F + GK p ) = 0 (300) Характеристическое уравнение наблюдателя состояния (299):

det[zE ( F K H C )] = 0, или после транспонирования выражения в квадратных скобках:

[ )] ( det zE F T C T K H = T Введем обозначения F Э = F T ;

G Э = GT ;

K p = K H.

Э T Тогда характеристическое уравнение наблюдателя принимает вид:

[ )] ( det zE F Э + G Э K p = Э (301) Сравнивая (300) и (301), убеждаемся, что задача синтеза наблюдателя состояния эквивалентна задаче синтеза эквивалентного П– регулятора:

r Эr u kЭ = K p xkЭ для эквивалентного объекта управления r r r x kЭ+1 = F Э x kЭ + Gu kЭ.

Решая эквивалентную задачу с помощью соотношений (289), (288), определяем матрицу коэффициентов передачи эквивалентного регулятора Э Kp и, следовательно, искомую матрицу коэффициентов передачи наблюдателя K H = K p.

Э Структурная схема многомерной дискретной АСР с ПИ – регулятором состояния – выхода и наблюдателем состояния приведена на рис. 81.

r Объект управления r r r uk yk x k +1 xk C Ee pT G F r yk Наблюдатель состояния KH r € r r yk € € xk +1 xk Ee pT G C Рис. F 81.

ПИ – регулятор Ee pT r € xk K r r y зад ek K 9. Многомерные дискретные АСР с прогнозом регулируемых переменных Перспективным способом решения дискретной задачи оптимального управления, рассмотренной в предыдущем разделе, является сведение её к последовательности задач математического программирования. Этот подход, в частности, позволяет:

– упростить вычислительную процедуру (вместо моделей объекта управления в пространстве состояния можно использовать более простые модели «вход – выход»;

не требуется решать уравнение Риккати);

– учитывать в явном виде ограничения на управления;

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

9.1. Структурная схема системы с прогнозом регулируемых переменных и его минимизацией [15] Структурная схема описываемой системы приведена на рис. 82.

Математическая модель многомерного объекта управления представлена в ней ( m n ) – матрицей переходных функций H (t ), связывающей n – r вектор регулируемых переменных y t с m – вектором кусочно–постоянных r управлений u t.

r r ut yt H (t ) r TУ rT € r yt +TУ yt И ut y зад Минимизация Прогноз Ф прогноза рассогласования r ut TУ Рис. 82.

Регулируемые переменные измеряются дискретно с интервалом дискретности TИ. Управления представляют собой кусочно–постоянные функции времени с периодом постоянства управления TУ. Предполагается, что TУ TИ.

Управляющее устройство системы включает два алгоритмических блока. В начале каждого шага управления первый блок осуществляет r прогноз вектора рассогласования y t между заданными и текущими значениями регулируемых переменных на величину текущего шага управления TУ. Второй блок осуществляет расчет изменения вектора r управлений на данном шаге ut, минимизирующего прогноз вектора r € рассогласования yt +TУ.

На изменение вектора управлений на каждом шаге могут быть наложены ограничения:

r r r u t u t u t+ (302) Итак, допустимое изменение вектора управлений на каждом шаге определяется в результате решения двух задач: прогноза рассогласования регулируемых переменных на текущий шаг управления и минимизации прогноза рассогласования в общем случае при наличии ограничений (302).

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

rr r ut = ut TУ + ut.

Найденное управление поддерживается неизменным до начала следующего шага управления фиксирующим элементом Ф.

9.2. Прогнозирование рассогласования [15, 16] r € Прогноз рассогласования yt +TУ может быть представлен суммой r € y tсвTУ, прогноза свободного (неуправляемого) рассогласования + вызванного действием на объект управления неконтролируемых возмущений, а также предыдущих управлений, и прогноза управляемого r упр € рассогласования yt +TУ, обусловленного ступенчатым изменением управления в начале текущего шага управления:

r r r € € € y t +TУ = y tсвTУ + y tупрУ (303) + +T Для прогноза свободного рассогласования регулируемых переменных применяются методы прогнозирования временных рядов. В частности, для получения краткосрочного прогноза в качестве прогнозируемой функции могут использоваться модели тренда временных рядов, для описания которого используют линейные или линеаризуемые относительно параметров функции. Например, в случае линейной прогнозирующей функции расчетные соотношения для прогноза, определяемые методом наименьших квадратов имеют вид:


r € y tсвTУ = b0 (t ) + b1 (t )TУ + 2 r n 1 n r (2n 1) yt i 3 i yt i, b0 (t ) = n(n + 1) i =0 i = 6 n 1 r r 2 n y t i n 1 i y t i, b1 (t ) = n(n + 1) i =0 i = где n – «память» прогнозатора – число членов временного ряда прогнозируемой переменной, используемых для прогнозирования.

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

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

Прогноз управляемого рассогласования определяется по выражению r r € y tупрУ = H (TУ )u t (304) +T 9.3. Минимизация прогноза рассогласования [15] В частном случае при отсутствии ограничений на управления (302) задача минимизации прогноза рассогласования может быть решена в явном виде. Максимальная точность регулирования достигается, если управления изменяются на каждом шаге так, чтобы обеспечить нулевое значение прогноза рассогласования:

r € yt +TУ = 0. (305) Подставляя в (305) выражения (303), (304) и разрешая полученное матричное уравнение, получаем выражения для оптимального управления при отсутствии ограничений:

[ ] r H (T )T H (T ) 1 H (T )T y св при m n € t +TУ r У У У 1 r св € ut = H (TУ ) yt +TУ при m = n [ ] T 1 r св € H (TУ ) H (TУ ) H (TУ ) yt +TУ при m n T mn Решения при (число управлений не равно числу регулируемых переменных) получены с использованием понятия обобщенной обратной матрицы и обладают минимальной нормой.

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

r yt +TУ min. (306) ut К сожалению, задача (306) с ограничениями (302) допускает только алгоритмическое решение. Одной из возможных постановок является решение задачи (306) как задачи с квадратичным критерием качества:

r r r r € € ytT+TУ Qyt +TУ + utT Rut min, (307) u t где Q и R – соответственно ( n n ) – и ( m m ) – диагональные весовые матрицы.

Первое слагаемое вводится в критерий задачи (307) для минимизации прогноза рассогласования, а второе – для минимизации изменения вектора управления.

Подставляя в (307) выражения (303), (304) и вводя обозначения:

r r r r r r v = ut u t ;

b = ut+ ut ;

D = H (Ty ) T QH (Ty ) + R, [ ] r r св r € c = 2 H (Ty ) Qyt +T + Du t T y можно преобразовать задачу (307) с ограничениями (302) к виду:

rr r r c Т v + v T Dv min r r (308) 0v b Задача (308) есть частный случай задачи квадратичного программирования, в которой требуется отыскать экстремум квадратичной целевой функции многих переменных при линейных ограничениях неравенствах и условиях неотрицательности переменных:

rr r r min Z = c T x + x T Dx, rrr (309) Ax b ;

x r где x - n-вектор переменных задачи;

r c - n-вектор коэффициентов линейных членов целевой функции Z;

D – (nn) – симметрическая матрица коэффициентов квадратичных членов целевой функции;

А – (mn) – матрица коэффициентов ограничений-неравенств;

r b - m-вектор свободных членов ограничений.

9.4. Сведение задачи квадратичного программирования к задаче о линейной дополнительности [17, 14] Одним из широко распространённых в настоящее время способов решения задачи квадратичного программирования является сведение этой задачи к задаче о линейной дополнительности.

Задачей о линейной дополнительности называется задача об r r отыскании векторов w и z таких, что выполняются условия:

r rr w = Mz + q (310.а) rT r w z = 0. (310.б) rr w, z 0 (310.в) Условия (310.а) есть ограничения-равенства. Входящие в них p rrr векторы w, z, q называются соответственно векторами базисных, свободных переменных и свободных членов ограничений.

M – (pp) – квадратная матрица коэффициентов ограничений.

Свободные переменные всегда равны нулю. Базисные переменные положительны. Значения базисных переменных при равенстве нулю свободных называют базисным решением системы (310.а).

Условие (310.б) называют условием дополняющей нежёсткости.

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

Условия (310.в) называются условиями неотрицательности переменных.

Для того чтобы свести задачу квадратичного программирования (309) к задаче о линейной дополнительности (310), запишем необходимые условия экстремума (условия Куна-Таккера) для задачи (309):

r rr r 2 = 2 Dx + AT 1 + c, rr r s = Ax + b, (311) rT r 2 x = 0,, rT r s 1 = 0, rrr r x, 1, 2, s 0 r r где 1, 2 - соответственно m- и n- векторы неопределённых множителей Лагранжа для ограничений-неравенств и условий неотрицательности переменных в задаче квадратичного программирования, r s - m-вектор дополнительных неотрицательных переменных, вводимых в ограничения-неравенства для преобразования их в равенства.

Вводя обозначения r r r 2 2 D | AT x c r r r M = | ;

w= ;

z = ;

q=, r r r 1 A | s b можем записать (311) в форме (310) и тем самым свести задачу квадратичного программирования к задаче о линейной дополнительности.

9.5. Решение задачи о линейной дополнительности методом Лемке r rr При q 0 задача (310) имеет тривиальное базисное решение: w = q.

r Если же вектор q содержит хотя бы один отрицательный элемент, начальное базисное решение оказывается недопустимым, т.к. нарушается условие (310.в). Для получения допустимого решения введём во все уравнения системы (310.а) искусственную положительную переменную z0, равную максимальному абсолютному значению среди всех отрицательных r элементов вектора q. Пусть для определённости min qi = qs, i тогда z 0 = qs.

При этом задача (310) принимает вид:

r rr r w M z e z 0 = q, (312.а) (312.б) wi zi = 0;

i = 1, p,, (312.в) wi, zi, z0 0 r где e - единичный вектор.

Для решения задачи (312) к системе ограничений-неравенств (312.а) применяют последовательность симплексных преобразований. Каждое симплексное преобразование позволяет поменять местами базисную и свободную переменные, т.е. базисную переменную вывести из базиса и превратить её в свободную, а свободную переменную ввести на её место в базис.

Симплексное преобразование легко формализовать при использовании симплексных таблиц. Начальная симплексная таблица для задачи (312) имеет вид (таблица 9) Таблица 9.

Базис w1 … ws … wl … wp z1 … zr … zp z0 q 1 -m1p - w1 -m11 -m1r q … ws -ms1 -msr -msp -1 qs … 1 -mlp - wl ql -ml1 -mlr … 1 -mpp - wp -mp1 -mpr qp В первом столбце таблицы записаны базисные переменные. Остальные столбцы соответствуют базисным, свободным переменным и искусственной переменной. Последний столбец содержит свободные члены ограничений-неравенств. На пересечении строк и столбцов, соответствующих базисным переменным, ставятся единицы, остальные элементы этих столбцов равны нулю. В столбцы, соответствующие свободным переменным, вписываются коэффициенты ограничений.

Элементы столбца, соответствующего искусственной переменной, полагаются равными минус единице.

Для того, чтобы свободную переменную zr ввести в базис, а базисную переменную wl вывести из базиса, необходимо осуществить симплексное преобразование относительно элемента mlr, находящегося на пересечении строки wl и столбца zr. Этот элемент и порождающие его строка и столбец называются разрешающими.

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

mir m mijП = mij mej ;

qiП = qi ir qe mer mer (индексом «П» снабжены значения коэффициентов и свободных членов ограничений после симплексного преобразования).

На первом шаге в базис вводится искусственная переменная z вместо базисной переменной ws (разрешающий элемент обозначен в таблице 9 кружком). После первого симплексного преобразования симплексная таблица принимает вид (таблица 10).

Таблица 10.

Базис w1 … ws … wp z1 … zs … zp z0 q 1 - w1 m1 p m1П q1П П m11 s … -1 z0 П msp msП П qsП mss … - wk П mkp mkП1 П qkП mks … -1 wp П П П П m p1 m ps m pp qp В результате начального симплексного преобразования базисная переменная ws превратилась в свободную, имеющую нулевое значение. В соответствии с условием дополняющей нежёсткости (312.б) дополнительная для переменной ws переменная zs должна иметь ненулевое значение, поэтому на втором шаге эту переменную следует ввести в базис.

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

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


qk q = min ( i, mis 0).

mks mis i На третьем шаге согласно правилу выбора разрешающего столбца в базис вводится переменная zk и т.д.

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

10. Автоматизация типовых технологических процессов [3, 18] 10.1. Регулирование основных параметров технологических процессов Регулирование расхода жидкости или газа Объектом регулирования является трубопровод. Скорость жидкости в трубопроводе определяется уравнением Бернулли:

p p v=c 2 = c 2g (313), где v – скорость жидкости, с1 – коэффициент расхода, g – ускорение силы тяжести,, =/g – соответственно удельный вес и плотность жидкости, p – перепад давления на трубопроводе.

Объёмный расход жидкости находим умножая скорость v на площадь сечения трубопровода f:

(314) Q=vf Уравнение статики трубопровода (баланс движущей силы потока Fдв и силы сопротивления трубопровода Fсопр):

Fдв = Fсопр, или с учётом (313), (314) Q fp = (315).

2 fc Из уравнения (315) можно найти коэффициент расхода Q c=.

f 2p Если приложенная к потоку сила Fдв превышает гидродинамическое сопротивление трубопровода, возникает ускорение потока dv/dt и, вместо уравнения статики, получаем уравнение dv + Fсопр = Fдв, m dt или с учётом (315), (314) m dQ Q + = fp, f dt 2 fc где m=V – масса жидкости в трубопроводе, V=Lf – объём трубопровода, L – длина трубопровода.

Подставляя в уравнение динамики выражение для массы, имеем:

dQ Q L + Q = fp.

dt 2 fc Наконец, приводя последнее уравнение к стандартному виду для инерционного звена первого порядка:

dQ + Q = kp T dt (T – постоянная времени, k – статический коэффициент передачи по каналу pQ), окончательно получаем:

2 Lfc 2 2 f 2c T= ;

k=.

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

Например, для L=60 м, диаметра трубопровода D=28 мм (f=0, м2), Q=40 л/мин (0,00066 м3/с), р=1,3 кгс/см2, =1000 кгс/м3, g=9,8 м/с постоянная времени трубопровода Т=0,5 с.

На практике находят применение три способа регулирования расхода.

1) Дросселирование потока на линии нагнетания (рис. 83) FC Рис. 83.

На рис. 83 обозначено:

- насос (компрессор), - рабочий орган с исполнительным механизмом, FC - регулятор (С) расхода (F).

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

2) Байпасирование – перепуск части потока из основного трубопровода в обводную линию (рис. 84).

FC Рис. 84.

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

3). Изменение напора в трубопроводе изменением числа оборотов вала насоса (рис. 85).

FC Рис. 85.

Здесь - регулируемый электропривод скорости вращения вала насоса.

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

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

Регулирование давления Объект регулирования – ёмкость с газом. Поведение идеальных газов описывается уравнением Менделеева-Клайперона:

p V = M R Ta, (316) где p – давление, V – объём, M = m µ - число молей газа, m – масса газа, µ - грамм/моль – число граммов равное молекулярному весу, R – универсальная газовая постоянная, равная работе, совершаемой единицей массы идеального газа при изобарном нагревании его на один градус, Ta – абсолютная температура.

Рассмотрим ёмкость с идеальным перемешиванием (Рис. 86):

G2, 2, p2= f2, p, G1, 1, p1 f1, Рис. 86.

Здесь G1, G2 – массовые расходы на линиях притока и расхода;

f1, f2 – проходные сечения клапанов, p1, p2, p – избыточные давления на входе, выходе и внутри ёмкости (считаем, что выходное давление p2 равно атмосферному);

1, 2, - плотности газа на входе, выходе и внутри ёмкости (в силу гипотезы об идеальном перемешивании 2 = ).

Уравнение статики ёмкости:

G1 = G2.

Уравнение динамики:

dm = G1 G2, dt или в приращениях:

dm = G1 G2. (317) dt Из уравнения (316) pVµ m=.

RTa При V, Ta = const имеем dm Vµ dp =. (318) dt RTa dt Согласно уравнениям (313), (314) объёмный расход газа равен p Q = cf 2, следовательно, массовый расход определяется выражением G = Q = cf 2 p.

В частности, расход на выходе ёмкости задаётся соотношением:

G2 = c2 f 2 2 p 2 (319) (поскольку p2 = 0).

Находя 2 = из уравнения Менделеева-Клайперона (316) m pµ 2 = =, V RTa подставляя последнее выражение в (319) и линеаризуя результирующее соотношение, находим µ G2 = c2 f 2 2 p. (320) RTa Подставляя (318) и (320) в (317), имеем Vµ dp 2µ + c2 f 2 p = G1, RTa dt RTa или в стандартной форме dp + p = kG1, T dt µ RTa V где T = k= ;

.

2µ c2 f c2 f 2 2 RTa Как видим, инерционность ёмкости пропорциональна её объёму, обратно пропорциональна проходному сечению f2 и корню квадратному из абсолютной температуры.

Дж кН Например, при V=50 м3, Та=20 оС, R=150, р=1000 2, f2=0, кг К м м2, С2=0,75 параметры модели объекта имеют значения: Т=115с, Н м k=10.

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

Регулирование уровня жидкости Объект регулирования – ёмкость с идеальным перемешиванием жидкости (рис. 87).

Q1, p1 f1, p f2, Q2, p2= Рис. 87.

Уравнение статики:

Q1 = Q2.

Уравнение динамики в приращениях:

dV (321) = Q1 Q2, dt где V = SH - объём жидкости.

При S = const dV dH =S (322).

dt dt Расход Q2 по-прежнему определяется уравнением Бернулли:

p p Q2 = c2 f 2 2, или при р2 = 0 и р=gH:

Q2 = c2 f 2 2 gH, следовательно, g Q2 = c2 f 2 H. (323) H (Индексом «0» помечены значения переменных, в окрестности которых осуществляется линеаризация нелинейных характеристик).

Подставляя (322), (323) в (321), имеем dH g + c2 f 2 H = Q1, S dt 2H или в стандартном виде:

dH + H = kQ1, T dt 2H 0 2H S где T = k= ;

.

c2 f 2 g c2 f 2 g Как видим, постоянная времени пропорциональна площади сечения, корню квадратному из уровня и обратно пропорциональна сечению сливного клапана.

Например, для ёмкости с характеристиками S=8 м2, f2=0,002 м2, м с2=0,6, Н0=2 м параметры модели равны: Т=4000с 1,1 час, k=500.

м3 с На практике находят применение следующие способы регулирования уровня.

1) Изменением расхода жидкости на входе в аппарат – регулирование на притоке (рис. 88) Qпр Н LC Рис. 88.

Здесь LC – регулятор уровня.

2) Изменением расхода на выходе аппарата – регулирование на стоке (рис. 89).

Н LC Qст Рис. 89.

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

3) Соотношением расходов на притоке и стоке (рис. 90) Н LC FFC Рис. 90.

Qст В данном случае для регулирования уровня используется каскадная АСР с промежуточной величиной – соотношением расходов на притоке и стоке (FFC – стабилизирующий регулятор соотношения расходов).

Каскадная АСР позволяет повысить качество регулирования уровня по сравнению с одноконтурными.

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

Например, АСР уровня жидкости в испарителе (рис. 91).

Регулирующее воздействие – подача теплоносителя в змеевик испарителя.

Пар Н LC Теплоноситель Рис. 91.

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

Экзотермический реактор с рубашкой (рис. 92) Gрс, рс р Gв1, в Gпр, пр Gв2, в Рис. 92.

Обозначения:

Gрс, Gпр – массовые расходы реакционной смеси и продукта, G в1, G в2 - массовые расходы охлаждающей воды на входе и на выходе рубашки, рс, пр, в1, в2,, р – температуры соответственно реакционной смеси, продукта, воды на входе и выходе рубашки, в реакторе и рубашке.

Допущения:

1) В реакторе и рубашке идеальное перемешивание:

пр = ;

в = р.

2) Потери тепла в окружающую среду отсутствуют.

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

4) Не учитывается динамика материальных потоков, т.е.

G р.с = Gпр ;

Gв = Gв = Gв.

1 5) Не учитывается тепловой эффект перемешивания.

6) Температура реакционной смеси и воды на входе в рубашку постоянны:

р.с = const ;

в = const.

С учётом допущения 3) объект является двухёмкостным.

Уравнение теплового баланса реактора:

d R рсV рс c рс = G рс с рс ( рс ) + 0 G рс s( р ). (324) µ dt Уравнение теплового баланса рубашки:

d р вVв cв = s ( р ) Gв св ( р в ), (325) dt рс, в – плотности реакционной смеси и воды, где Vрс,Vв – объёмы реакционной зоны реактора и рубашки, срс, св – удельные теплоёмкости реакционной смеси и воды, R0 – количество тепла выделяемое в процессе реакции одного моля реакционной смеси, µ - молярный вес реакционной смеси, - коэффициент теплопередачи через стенку между реактором и рубашкой, s – поверхность теплопередачи.

Записывая уравнения (324), (325) в приращениях и приводя их к стандартному виду, получаем систему двух уравнений первого порядка, описывающих рассматриваемый объект:

d + = K1G рс + K 2 р T dt, (326) d р + р = K 3 + K у Gв T dt где постоянные времени и коэффициенты передачи определяются соотношениями:

рсV рс с рс вVв св T1 = T2 = ;

;

G рс, 0 с рс + s Gв, 0 св + s R с рс ( рс 0 ) + s µ K1 = K2 = ;

;

G рс, 0 с рс + s G рс, 0 с рс + s c в ( р, 0 в ) s K4 = K3 = ;

.

Gв, 0 с в + s Gв, 0 с в + s Передаточная функция объекта по каналу Gв р :

K W1 ( p) =.

T2 p + Передаточная функция объекта по каналу р :

K W2 ( p) =.

T1 p + Основной задачей регулирования теплового режима реактора является поддержание температуры в реакционной зоне реактора изменением расхода воды в рубашку. Передаточная функция объекта по каналу регулирующего воздействия Gв :

K2K Wx ( p ) = W1 ( p )W2 ( p ) =.

(T1 p + 1)(T2 p + 1) Передаточная функция объекта по каналу возмущающего воздействия (колебаний расхода реакционной смеси) G рc :

K WF ( p ) =.

T1 p + Теплообменник типа «труба в трубе»

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

Теплообменник представляет две концентрические трубы (рис. 93).

Нагреваемая жидкость подаётся во внутреннюю трубу, теплоноситель – во внешнюю.

Обозначения:

x – координата (длина трубы), Dв – диаметр внутренней части внутренней трубы, - толщина стенки внутренней трубы, Dв + 2 - диаметр внешней части внутренней трубы, Dн – внутренний диаметр наружной трубы, L1 = Dв – длина внешней окружности внутренней трубы, L2 = (Dв + 2) – длина внешней окружности внутренней трубы, L3 = Dн – длина внутренней окружности наружной трубы, Sв - Dв2/4 – сечение внутренней трубы, Sст = [(Dв + 2)2-Dв2]/4 = ( Dв + ) – сечение стенки внутренней трубы, Sк = [Dн2 - (Dв + 2)2]/4 – сечение кольца между внутренней и внешней трубами.

теплоноситель нагреваемая Dв Dн жидкость Рис.93.

x x Объект содержит три тепловые ёмкости: тепловая ёмкость нагреваемой жидкости во внутренней трубе, ёмкость стенки внутренней трубы, ёмкость нагревателя в межтрубном пространстве.

Для получения уравнения динамики изменения температуры материала во внутренней трубе рассмотрим элемент объёма жидкости внутренней трубы Sв·х (х – приращение длины трубы относительно значения х).

Количество тепла в элементарном объёме Sвх:

Sвхжсжж, где ж, сж, ж – соответственно плотность, удельная теплоёмкость и температура нагреваемой жидкости.

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

Приращение тепла в элементе Sвх за промежуток t:

ж S в x ж сж t. (327) t Найдём приращение тепла в элементе объёма за счёт приращения координаты х. Поток тепла через сечение внутренней трубы Sв за счёт перемещения жидкости:

жсжжvжSв, где vж – скорость движения жидкости.

Количество тепла через сечение Sв за промежуток t:

жсжжvжSвt.

Приращение тепла по длине теплообменника:

ж ж сж vж S в tx. (328) x Сумма (327) и (328) есть общее приращение тепла в элементе объёма Sвх за промежуток времени t.

Тепло, отдаваемое стенкой трубы элементу Sвх за время t:

1L1(ст-ж)xt, 1 – коэффициент теплоотдачи от стенки трубы нагреваемому где материалу, L1x – поверхность теплоотдачи, cт – температура стенки.

Уравнение теплового баланса для нагреваемой жидкости:

ж ж S в ж сж + ж сж S в vж = 1 L1 ( ст ж ) (329) t x Разделим левую и правую части уравнения (329) на Sвжсж и запишем его в виде:

ж ж + vж + c1 ж = c1 ст, (330) t x где 1 L c1 =.

S в ж сж Динамика изменения температуры теплоносителя во внешней трубе описывается уравнением аналогичным (330):

т т + vт + c2 т = c2 ст, (331) t x где т, ст, т – параметры теплоносителя, vт – скорость теплоносителя, 2 L c2 =, S к т ст 2 – коэффициент теплопередачи от теплоносителя к стенке.

Для описания динамики изменения температуры стенки рассмотрим элементарный объём стенки Sстх. Изменение количества тепла, запасённого элементом объёма стенки:

ст S ст x ст, t где ст, сст, ст – параметры стенки.

Уравнение теплового баланса стенки:

ст S ст cст ст = 2 L2 ( т ст ) 1 L1 ( ст ж ), t или ст + (с3 + с4 ) ст = с3 т с4 ж, (332) t где 2 L2 1 L c3 = c4 = ;

.

S ст ст сст S ст ст сст Уравнения (330), (331), (332) образуют систему трёх уравнений, описывающих динамику теплообменника «труба в трубе». Для получения зависимости температуры нагреваемой жидкости от температуры теплоносителя необходимо решить эту систему уравнений при заданных начальных условиях.

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

10.2. Типовые схемы автоматизации технологических процессов Автоматизация насосов и компрессоров Насосы и компрессоры предназначены для перемещения соответственно жидкостей и газов.

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

изменение вязкости и плотности перемещаемой жидкости.

Схема автоматизации насоса (компрессора) показана на рис. 94.

Обозначения:

I, R, C – показывающий, регистрирующий, регулирующий приборы;

А – световая и звуковая сигнализация;

Н – ручное управление;

J – многоточечный (обегающий) прибор;

P, F (в первой позиции), T, E – соответственно давление, расход, температура, электрические показатели.

Как отмечалось в разделе 10.1, наиболее простым способом регулирования расхода является дросселирование потока жидкости (газа) на линии нагнетания (позиция 1). Схема автоматизации также предусматривает: контроль давления на всасывающей и нагнетающей линиях (поз. 2, 3) и сигнализацию на нагнетающей линии;

контроль и сигнализацию температуры хладагента, подшипников и обмоток электродвигателя насоса (поз. 4);

контроль и сигнализацию давления хладагента и смазки (поз. 5, 6) и расхода хладагента и смазки (поз. 7, 8);

сигнализацию отклонения режима работы двигателя от нормального – перегрузку двигателя (поз. 9). Предусмотрены ручное управление задвижками на линиях всасывания и нагнетания (поз. 10, 11) и сигнализация положения задвижек (поз. 12, 13).

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

Автоматизация поверхностных кожухо-трубчатых теплообменников (рис. 95) Предназначены для нагрева продукта теплом от теплоносителя.

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

Основной показатель качества – температура нагреваемого продукта – стабилизируется расходом теплоносителя (поз.1). Основные возмущения:

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

Теплообменники (и тепловые объекты вообще) характеризуются значительной инерционностью.

Контролируются: расходы теплоносителя, (поз.2) и продукта (поз.3), их конечные и начальные температуры (поз.4), давление теплоносителя и продукта (поз.5,6).

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

Рис. 95.

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

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

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



Pages:     | 1 | 2 || 4 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.