Данная страница содержит автореферат моей будущей магистерской работы. Вполне возможно, что представленная к защите работа будет сильно отличаться от ниже следующего.
Виталий Орехов
АВТОМАТИЗАЦИЯ СИНТЕЗА РЕГУЛЯТОРОВ И НАБЛЮДАТЕЛЕЙ СОСТОЯНИЯ
В СРЕДЕ ПАКЕТА MATLAB
Синтез регуляторов состояния (РС) и наблюдателей состояния (НС) является
весьма трудоемкой задачей, решение которой можно упростить, используя современное
программное обеспечение, предназначенное для автоматизации математических и
научных расчетов, например, пакеты Mathematica, MathCAD, MatLab, Maple. Система
программирования MatLab выгодно отличается от других "математических" пакетов
наличием в ней приложения для структурного математического моделирования Simulink,
инструментов, специально предназначенных для анализа и синтеза линейных систем
автоматического управления (Control System Toolbox) и развитого объектно-
ориентированного алгоритмического языка, а также возможностью доступа к своим
функциям. Однако, чем больше возможностей предоставляет пользователю программный
пакет, тем более необходимыми становятся методические рекомендации по их
рациональному использованию при решении конкретных инженерных и научных задач.
1 МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ
Исходная (а) и преобразованная (б) матричные структурные схемы системы
автоматического управления (САУ) с регулятором состояния полного порядка (РСПП)
представлены на рис. 1.1, а структурные схемы наблюдателя состояния полного
порядка (НСПП) - на рис. 1.2. Сразу оговоримся, что почти все описываемые в
данной работе процедуры синтеза справедливы только для систем с единственным
управляющим сигналом (т.е. для одномерных систем).
а) б)
Рис. 1.1. Структурные схемы САУ с регулятором состояния полного порядка
а) б)
Рис. 1.2. Структурные схемы наблюдателя состояния полного порядка
Синтез РСПП и НСПП состоит в выборе коэффициентов корректирующих обратных
связей K и L соответственно из условия обеспечения желаемого распределения
полюсов или, что одно и то же, желаемых коэффициентов характеристического
полинома замкнутой САУ с РСПП или замкнутого НСПП.
Таким образом, для определения искомых коэффициентов необходимо решить одно
из уравнений
(1.1)
или
, (1.2)
где Aз - матрица состояния замкнутой системы, которая определяется для САУ
с РСПП как
для НСПП как
- вектор желаемых полюсов замкнутой системы;
- желаемый характеристический полином замкнутой системы.
Существуют также аналитические выражения для вычисления искомых матриц
коэффициентов РСПП и НСПП. Так, в соответствии с опубликованной в 1980г. формулой
Аккермана [3] коэффициенты РСПП можно определить из последней строки матрицы
(1.3)
где
- матрица управляемости.
Ещё один способ основывается на использовании прямого и обратного
преобразования математического описания разомкнутой и замкнутой систем в
пространстве состояний к канонической форме Фробениуса. Согласно алгоритму
синтеза, приведенному в [1, 2], уравнения
и
имеют аналитические решения:
(1.4)
(1.5)
где
- коэффициенты характеристического полинома разомкнутой системы (или разомкнутого
наблюдателя);
- коэффициенты желаемого характеристического полинома;
- матрица наблюдаемости;
, .
Другой способ поиска вектора искомых коэффициентов требует представления
математического описания системы управления в канонической форме фазовой
переменной, или перехода к описанию в канонической форме, если исходные уравнения
состояния записаны в естественном базисе [1]. Пара матриц А и B считается
записанной в каноническом базисе, если эти матрицы имеют вид:
,
Следует отметить, что пара матриц А и А* в естественном и каноническом базисах
имеют одни и те же собственные значения, кроме того, существует матрица
преобразования (перехода) к канонической форме Р, такая, что выполняются условия:
,
Неособая матрица Р при этом задаётся в виде:
где
Представив описание объекта в канонической форме фазовой переменной,
вычисляют вектор коэффициентов регулятора состояния в каноническом базисе как
(1.6)
где - коэффициенты характеристических полиномов соответственно объекта
(разомкнутой системы) и желаемого (замкнутой системы).
Имея матрицу перехода к каноническому базису Р, можно получить действительные
коэффициенты обратных связей
(1.7)
Методики, описываемые формулами (1.6, 1.7) и (1.4, 1.5), удобны тем, что
в выражения входят разности коэффициентов при соответствующих степенях
характеристических полиномов разомкнутой системы и желаемого стандартного
полинома. А все остальные задачи - формирование матриц и матричные операции,
предусмотренные тем или иным методом синтеза - хорошо поддаются алгоритмизации,
и не представляют затруднений.
Формулы (1.1) и (1.2), (1.3), (1.4) и (1.5), (1.6) и (1.7) могут с успехом
лечь в основу как численных, так и аналитических методов автоматизации синтеза
регуляторов и наблюдателей состояния.
Используемые в аналитических методах чинтеза РСПП и НСПП матрицы
преобразования при высоком порядке дифференциальных уравнений, описывающих объект
регулирования, могут оказаться плохо предопределенными и давать большую
погрешность. В связи с этим заслуживают внимания также и численные алгоритмы
синтеза регуляторов и наблюдателей состояния полного порядка, наиболее популярный
из которых предложен в [4].
Следует отметить, что, если известен какой-либо алгоритм поиска коэффициентов
РСПП, то подставляя в него вместо матриц A и B матрицы AТ и CТ соответственно,
можно определить транспонированный вектор коэффициентов НСПП LТ. Эта особенность
обусловлена дуальной связью между управляемостью и наблюдаемостью.
Итак, при описании алгоритмов синтеза РСПП и НСПП мы будем рассматривать
четыре основных способа:
1. Алгоритм, основанный на составлении и решении системы n линейных
алгебраических уравнений с n неизвестными относительно искомых коэффициентов
(в матричной трактовке этот алгоритм описывается формулой (1.2)).
2. Алгоритмы расчета, использующие стандартные (т.е. входящие в стандартную
поставку системы MatLab) функции. Это прежде всего, способ, основанный на
использовании формулы Аккермана (1.3), и численный метод, предложенный в [4].
3. Алгоритм, описанный в [1,2] и основанный на использовании формул (1.4) и (1.5).
4. Алгоритм вычислений с использованием перехода к описанию объекта в
канонической форме фазовой переменной (формулы (1.6),(1.7)).
2 МЕТОДИКА МАШИННОГО СИНТЕЗА РЕГУЛЯТОРОВ
СОСТОЯНИЯ И НАБЛЮДАТЕЛЕЙ СОСТОЯНИЯ ПОЛНОГО
ПОРЯДКА В ЧИСЛЕННОМ ВИДЕ
Первый способ. В численном виде этот метод плохо поддается автоматизации,
поскольку требует от разработчика усилий по составлению системы алгебраических
уравнений. При составлении системы правильность решения задачи и затраченное время
зависит от умения исследователя быстро и безошибочно делать структурные
преобразования, выводить передаточные функции и упрощать громоздкие аналитические
выражения. Полученную систему уравнений затем можно решить одним из предлагаемых
системой MatLab способов. Между тем, задача составления системы уравнений намного
упрощается при использовании средств символьной математики. Поэтому мы вернёмся к
этому способу в следующем разделе.
Второй способ. Система программирования MatLab предлагает пользователю две
функции для расчета коэффициентов регулятора состояния по желаемым полюсам [5]:
K=acker(A,B,Р)
и
K=place(A,B,Р)
Обе решают одну и ту же задачу (формула (1.1)) разными методами: первая - по
формуле Аккермана (1.3), а вторая - по численному алгоритму, предложенному в [4].
Разработчики MatLab рекомендуют применять функцию acker для синтеза систем с одним
входом и одним выходом, описываемых дифференциальными уравнениями не выше 5-го
порядка, а функцию place - во всех остальных случаях, включая и многомерные
системы, хотя нам удавалось во многих случаях успешно применять функцию acker и
для систем более высокого порядка. Функция acker по окончании работы проверяет
результат синтеза, выполняя попарное сравнение значений желаемых и полученных в
результате синтеза полюсов системы. Перед сравнением полюсы сортируются. Следует
отметить небольшую неточность, допущенную авторами пакета при разработке этой
функции. Применяемый в функции acker способ сортировки (sort) дает неоднозначный
результат в случае комплексно-сопряженных чисел с одинаковой амплитудой, что
имеет место, например, в случае использования полинома Баттерворта. В результате,
несмотря на абсолютно точно найденный результат, на экран выводится предупреждение
"Warning: Pole locations are more than 10% in error". Для ликвидации этого недоразумения
необходимо заменить в функции acker метод sort методом esort. Исправленную
функцию можно записать под другим именем, например, acker_m,поскольку в новых
версиях системы MatLab изменения существующих функций не ограничиваются просто
изменениями в текстовых m-файлах. Листинг изменённой функции приведен ниже.
function k = acker_m(a,b,p)
error(nargchk(3,3,nargin));
[msg,a,b]=abcdchk(a,b); error(msg);
[m,n] = size(b);
if n ~= 1 error('System must be single input');
end
p = p(:); % Make sure roots are in a column vector
[mp,np] = size(p);
[m,n] = size(a);
if (m ~= mp) error('Vector P must have SIZE(A) elements');
end
% Form gains using Ackerman's formula
k = ctrb(a,b)\polyvalm(real(poly(p)),a);
k = k(n,:);
% Проверка правильности численного результата
% (методы сортировки sort заменены на esort)
p = esort(p);
i = find(p ~= 0);
p = p(i);
pc = esort(eig(a-b*k)); pc = pc(i);
if max(abs(p-pc)./abs(p)) > .1
disp('Warning: Pole locations are more than 10% in error.')
end
Третий способ. Для решения поставленной задачи в соответствии с данным
методом можно предложить использовать дополнительно написанную функцию - например,
такую, текст которой приведен ниже. Функция kuo организует вычисление коэфициентов
регулятора состояния по формуле (1.4). Следует заметить, что в отличие от acker и
place, работающих с желаемыми полюсами, функция kuo оперирует коэффициентами
желаемого полинома.
function K=kuo(A,B,g)
n=length(g)-1;
a=real(poly(A)); %Команда получения коэфф-ов полинома объекта
M=eye(n,n);
for i=2:n
for j=1:i-1
M(i,j)=a(i-j+1);
end
end
Qy=ctrb(A,B);
Ks=g-a;
Ks=Ks(2:n+1);
K=((M*Qy')\Ks')';
Четвёртый способ. Для реализации вычислений в соответствии с формулами (1.6)
и (1.7) необходимо написать функцию, пример которой приведен ниже. Данная функция
также, как и функция kuo, использует коэффициенты желаемого полинома.
function K=trans_rf(A,B,g)
Qy=ctrb(A,B);
n=length(Qy);
P1=zeros(1,n);
P1(n)=1;
P1=P1*Qy^(-1);
P(1,:)=P1;
for i=2:n
P(i,:)=[P1*A^(i-1)];
end
a=real(poly(A)); %Команда получения коэфф-ов полинома объекта
K_canon=g-a;
K_canon=K_canon(2:(n+1));
K=fliplr(K_canon)*P;
Комментарии, касающиеся всех описанных способов. Как уже было сказано выше,
вследствие дуальной связи между наблюдаемостью и управляемостью, любую из
приведенных функций можно использовать и для синтеза НСПП, например,
Lt=acker_m(A', C', P); L=Lt'
(операция A' в пакете MatLab означает транспонирование матрицы).
Система MatLab содержит также функции вычисления полюсов некоторых
стандартных полиномов, например, Баттерворта, Бесселя, Чебышева [6-8]. Они
содержатся в папке MatLab\toolbox\signal (Signal Processing Toolbox) [9].
Определить полюсы нормированного по среднегеометрической частоте полинома
Баттерворта n-го порядка можно при помощи функции
[z,P,k]=buttap(n),
а ненормированного (с заданным среднегеометрическим корнем w) - при помощи функции
[z,P,k]=butter(n,w,'s').
По найденным полюсам можно определить и коэффициенты полинома gamma,
используя команду
gamma=real(poly(P)).
Для конструирования других стандартных полиномов пользователь может написать
собственные функции.
Например, функция для определения полюсов и коэффициентов биномиального
полинома [7], может иметь следующий вид:
function [P,g]=bin(n,w)
P=w*ones(1,n);
g=real(poly(P));
Одним из достоинств пакета MatLab я вляется то, что он может помочь
исследователю определить математическое описание объекта регулирования в
пространстве состояний в численном виде по его структурной схеме. Чтобы
воспользоваться этой возможностью, необходимо в приложении Simulink набрать
структурную модель объекта, пометив ее вход и выход портами, записать эту модель в
файл (*.mdl) и выполнить в программном режиме или из командной строки MatLab
команду
[A,B,C,D]=linmod('имя-файла')
Для конкретизации места расположения переменных состояния в структурной
схеме Simulink-модель желательно набрать в детализированной форме. Например, для
объекта управления "Тиристорный преобразователь - двигатель постоянного тока",
имеющего известную структурную схему, приведенную на рис. 2.1,а, Simulink-модель,
необходимая для использования функции linmod при синтезе регулятора состояния,
может иметь вид рис. 2.1,б.
а)
б)
Рис. 2.1 Структурные схемы объекта "Тиристорный преобразователь - двигатель
постоянного тока".
При реальном замыкании обратных связей после их расчета, следует определить
очередность следования переменных состояния при помощи команды
[[sizes,x0,xstr]=имя_файла_модели
Интересующая нас информация будет содержаться в переменной xstr в виде путей
к переменным состояния, например,
xstr =
'TPD/Integrator1'
'TPD/Integrator2'
'TPD/Integrator'
3 МЕТОДИКА МАШИННОГО СИНТЕЗА РЕГУЛЯТОРОВ
СОСТОЯНИЯ И НАБЛЮДАТЕЛЕЙ СОСТОЯНИЯ ПОЛНОГО
ПОРЯДКА В СИМВОЛЬНОМ ВИДЕ
Очень часто исследователя интересуют не только численные значения
коэффициентов регулятора и наблюдателя состояния, но и их аналитические выражения.
Современные математические пакеты могут справиться и с этой (гораздо более сложной
для ЭВМ) задачей за счет включения в них инструментов символьной математики.
В системе MatLab эти функции выполняет Extended Symbolic Toolbox [7], основу
которого составляет ядро системы символьной математики Maple.
Первый способ. Решение системы n уравнений с n неизвестными. При решении
поставленной задачи "вручную" обычно поступают следующим образом:
- находят выражения для характеристического полинома замкнутой САУ или замкнутого
НС через параметры разомкнутого объекта регулирования и коэффициенты обратных
связей;
- составляют систему n линейных уравнений с n неизвестными, приравнивая
коэффициенты при одинаковых степенях оператора Лапласа действительного и
желаемого характеристического полиномов;
- решают полученную систему относительно коэффициентов ki или li регулятора или
наблюдателя состояния.
Провести машинный синтез в соответствии с этим методом можно только в
символьном виде. Более того, исследователи наталкиваются на сложность решения
этой задачи в такой трактовке и вынуждены использовать для её решения другие
прикладные пакеты,например, Maple. Покажем способ решения с использованием системы
MatLab на простом примере. Синтезируем РСПП для системы "Тиристорный
преобразователь - двигатель постоянного тока", представленной на рис. 2.1.
Объявим все необходимые символьные переменные, используя описатель syms:
syms Tmu Tm C Ra Kd Kp Ta K k1 k2 k3 A B p w0
Зададим описание данного объекта в виде матриц пространства состояний. Эти
матрицы выводятся вручную, в соответствии с системой дифференциальных уравнений,
или по структурной схеме.
A=[0, Ra/C/Tm, 0; -C/Ra/Ta, -1/Ta, 1/Ra/Ta; 0, 0, -1/Tmu];
B=[0; 0; Kp/Tmu];
Получить характеристический полином замкнутой через РСПП системы можно по
известному выражению:
Hn=det(p*I-A+B*K);
где I - диагональная единичная матрица (получить её можно командой eye),
или, что ещё проще
syms p, Hn=poly(A-B*K,p)
При необходимости полученное выражение можно упростить с помощью таких
функций, как simple, collect, simplify. Для более наглядного вывода на экран
результатов можно пользоваться функцией pretty.
Для описываемого примера результат будет иметь вид:
2
3 (Tm Tmu + Tm Ta) p (Tmu + Tm) p 1
p + ------------------- + -------------- + -----------
Ta Tm Tmu Ta Tm Tmu Ta Tm Tmu
Далее необходимо из полученного символьного объекта - характеристического
полинома сформировать вектор его коэффициентов. К сожалению, система MatLab не
обладает готовым решением этой задачи, т.е. нет стандартной функции, работающей с
символьными полиномами. Имеющиеся в системе MatLab функции poly2sym и sym2poly
работают только с полиномами, у которых коэффициенты заданы в численном виде.
Функцию преобразования полинома в вектор символьных коэффициентов несложно
создать самостоятельно. Предлагаемый вариант такой функции приведен ниже. При
вызове данной функции ей необходимо передать как параметры полином и переменную,
по степеням которой данный полином составлен (например, в случае непрерывных
систем - это обычно переменные p или s, в случае дискретных систем - z).
function [a]=poly2koef_m(H,p)
i=1;
while H ~= 0
a(i)=subs(H,p,0);
temp=H-a(i);
H=collect(temp, p);
H=simple(H/p);
i=i+1;
end
a=fliplr(a);
При наличии такой функции коэффициенты характеристического полинома
замкнутой системы в виде вектора символьных объектов можно получить командой:
a_nz=poly2koef_m(Hn,p);
Коэффициенты полинома, соответствующего желаемому распределению корней,
также задаются в виде вектора символьных объектов или в виде вектора числовых
значений. Например, для желаемого полинома третьего порядка такой вектор в
символьном виде может иметь вид:
g=[1 g1*w0 g2*w0^2 g3*w0^3];
Теперь нужно составить систему уравнений, приравнивая коэффициенты
действительного и желаемого полиномов, и решить её относительно вектора
коэффициентов РСПП. Сделать это можно с помощью следующей последовательности
команд:
i=[1:3]; % Система 3 - го порядка
eqn(i)=a_nz(i+1)-g(i+1); % Составляем систему уравнений
[k1 k2 k3]=solve(eqn(1),eqn(2),eqn(3),k1,k2,k3)
K=[k1 k2 k3].'
Результат будет следующим (для вывода пользуемся функцией pretty(K)):
[ 3 2 ]
[ C Tmu (-1 + g1 w0 Ta - g3 w0 Tm Ta ) ]
[ - --------------------------------------- ]
[ Kn Ta ]
[ 2 2 ]
[ Tmu Ra (Ta - Tm + Tm g1 w0 Ta - g2 w0 Ta Tm) ]
[- ------------------------------------------------ ]
[ Kn Tm Ta ]
[ ]
[ -Tmu - Ta + g1 w0 Ta Tmu ]
[ ----------------------------- ]
[ Kn Ta ]
Второй способ. Для того, чтобы воспользоваться формулой Аккермана (1.3),
нужно "заставить" функцию acker работать с символьными объектами. Для этого
необходимо исключить из неё ту часть, которая отвечает за проверку правильности
произведенных вычислений. Поскольку собственно вычислительная часть функции
использует простые матричные операции, то полученная функция сможет успешно
работать и с символьными объектами тоже. Кроме того, при решении в символьном виде
гораздо удобнее задавать не желаемые полюсы p, а коэффициенты желаемого полинома g,
представленные в виде символьного вектора. В этом случае придётся переделать
функцию ещё раз, заменив строку
k = ctrb(a,b)\polyvalm(real(poly(p)),a);
где используется команда получения коэффициентов полинома real(poly(p)) на строку
k = ctrb(a,b)\polyvalm(g,a);
а в строке объявления функции заменить параметр р (желаемые полюсы) параметром g
(коэффициенты желаемого полинома). Полученную функцию необходимо сохранить под
каким-либо другим именем, например, acker_sym.
Тогда для рассматриваемого примера вектор искомых коэффициентов можно найти
с помощью команды:
K=acker_sym(A,B,g);
Результат в данном случае получается аналогичный (см. выше).
Третий и четвёртый способы. Для работы в символьном виде функции kuo и trans_rf,
описанные в предыдущем разделе, также нуждаются в некоторой доработке. В этих
функциях необходимо команду получения коэффициентов полинома объекта (листинг
этих функций см. выше):
a=real(poly(A));
заменить на следующую последовательность команд:
syms p % оператор Лапласа;
H=poly(A,p); % формирование полинома;
H=collect(simple(H),p); % упрощение;
a=poly2koef_m(H,p); % формирование вектора коэффициентов,
Если изменённые функции сохранить под именами kuo_sym и trans_rf_sym, то
с помошью команд
K=kuo_sym(A, B, g);
или
K=trans_rf_sym(A, B, g);
мы получим уже знакомый результат (см. выше).
4 ФОРМИРОВАНИЕ СТАНДАРТНЫХ ДИСКРЕТНЫХ
ПОЛИНОМОВ ПРИ СИНТЕЗЕ РЕГУЛЯТОРОВ
И НАБЛЮДАТЕЛЕЙ СОСТОЯНИЯ
Первым этапом синтеза цифровых регуляторов и наблюдателей состояния является
выбор желаемого стандартного характеристического полинома и преобразование его в
дискретную форму. Коэффициенты дискретного полинома зависят от периода прерывания
и типа корней соответствующего аналогового прототипа. С повышением порядка систем
сложность и количество соотношений, связывающих коэффициенты дискретного полинома
с корнями аналогового полинома, возрастает. Для упрощения и ускорения этого этапа
синтеза целесообразна его автоматизация.
Пусть известны коэффициенты некоторого аналогового полинома n-ой степени
, нормированного по его среднегеометрическому корню (СГК):
, (4.1)
, (4.2)
, (4.3)
. (4.4)
Один из наиболее удобных способов дискретизации полинома (4.1) состоит из
следующих этапов:
1) определение корней полинома (4.3) и корней полинома (1):
; (4.5)
2) определение корней дискретного полинома:
, (4.6)
где Т - период дискретности;
3) формирование дискретного полинома
. (4.7)
В пакете MatLab информация о полиномах хранится в виде массива
коэффициентов, упорядоченных по убыванию степени независимой переменной. Например,
вектор
содержит информацию о полиноме (4.3).
Корни полинома по его коэффициентам определяет функция roots, например,
P=roots(G), а обратную операцию выполняет функция poly, например, G=real(poly(P)).
Следует отметить, что в справочной литературе могут быть даны как
коэффициенты, так и корни стандартных полиномов, однако чаще приводятся значения
коэффициентов . Поэтому сведем в табл. 1 значения нормализованных корней
стандартных полиномов, наиболее часто используемых при синтезе систем управления
электромеханическими объектами: 1 - полином Баттерворта [6-8], 2 - полином
с 5%-ным перерегулированием [11], 3 - полином Грехема-Летропа [12, 8], 4 - полином
Бесселя [7], 5 - полином, сформированный методом двойных пропорций [13],
6 - полином с кратными комплексно-сопряженными корнями и коэффициентом
демпфирования 0,75. В таблице отсутствуют корни полиномов с биномиальными
коэффициентами, поскольку все они равны -1.
Таблица 1 - Нормализованные корни некоторых стандартных полиномов
Поли- ном
1
2
3
4
5
6
n=2
-0,707±0,707i
-0,689±0,724i
-0,700±0,714i
-0,866±0,500i
-0,707±0,707i
-0,750±0,661i
n=3
-0,500±0,866i -1
-0,571±0,821i -1
-0,521±1,068i -0,708
-0,746±0,711i -0,942
-0,500±0,866i -1
-0,750±0,661i -1
n=4
-0,383±0,924i -0,924±0,383i
-0,501±0,865i -0,940±0,342i
-0,424±1,263i -0,626±0,414i
-0,657±0,830i -0,905±0,271i
-0,707±0,707i -0,707±0,707i
-0,750±0,661i -0,750±0,661i
n=5
-0,309±0,951i -0,809±0,588i -1
-0,456±0,890i -0,853±0,522i -1
-0,376±1,292i -0,376±1,292i -0,896
-0,591±0,907i -0,852±0,443i -0,926
-0,378±0,441i -1,122±1,307i -1
-0,750±0,661i -0,750±0,661i -1
Для формирования стандартных полиномов, вычисление параметров которых не
предусмотрено в базовом наборе функций системы MatLab, можно написать собственную
функцию, восполняющую этот пробел. Например, корни и коэффициенты полинома с
биномиальными коэффициентами сможет определять такая MatLab-функция пользователя:
function [P ,a] = bin (n, w0)
if nargin == 1, w0 = 1; end
P = w0 * ones (1, n);
a = real (poly (P) );
Операции, заданные формулами (4.5), (4.6), в среде MatLab нет необходимости
выполнять в цикле, так как в алгоритмическом языке этого пакета предусмотрено
выполнение матричных операций, например,
z = exp (T * P).
Более того, операции над матрицами в целом здесь выполняются быстрее, чем
последовательность операций над их элементами.
Таким образом, фрагмент MatLab-программы, вычисляющей коэффициенты
дискретного полинома, соответствующего аналоговому полиному Бесселя 3-го порядка,
может выглядеть следующим образом:
w0 = input ('Введите СГК ')
[b,a] = besself (3, w0)
P = roots (a)
T=1
Z = exp (T * P)
ad = poly (Z)
Листинг выполнения этой программы имеет вид
Введите СГК 1
w0 =
1
b =
0 0 0 1
a =
1.0000 2.4329 2.4662 1.0000
P =
-0.7456 + 0.7114i
-0.7456 - 0.7114i
-0.9416
T =
1
Z =
0.3594 + 0.3097i
0.3594 - 0.3097i
0.3900
ad =
1.0000 -1.1087 0.5054 -0.0878
откуда следует, что аналоговый полином Бесселя 3-го порядка с единичным СГК
и соответствующий ему дискретный полином при определяются следующими
выражениями:
,
.
Несколько сложнее поставленная задача решается в символьном виде. В этом
случае нам необходимо найти не численные значения коэффициентов дискретного
полинома, а аналитические зависимости этих коэффициентов от параметров корней его
аналогового прототипа. Решение задачи усложняется тем, что вид указанных
аналитических выражений зависит от того, сколько корней образуют комплексно-
сопряженные пары и сколько из них являются действительными, а также от
соотношения амплитуд и фаз, действительных и мнимых составляющих этих корней.
Вариант программы, решающей поставленную задачу для полиномов 2-го порядка,
может иметь такой вид:
syms z p a b T
P(1)=a+i*b; P(2)=w*a-i*b
Z=exp(T*P);
G=prod(z-Z);
G=expand(G);
G=subs(G,exp(i*T*b),cos(T*b)+i*sin(T*b));
G=simplify(G);
G=collect(G,z)
G_bin=subs(G,[a b],[-w 0])
T=0.005, w=1/T
G_bin_c=sym2poly(eval(G_bin))
Описатель syms объявляет следующие за ним переменные символьными.
Для работы с комплексными числами в пакете MatLab предопределены
переменные i и j, значением которых является . Таким образом, в программе
первоначально заданы два комплексно-сопряженных корня
. (4.8)
Функция prod вычисляет произведение всех элементов вектора-аргумента. Из
функций преобразования символьных выражений в данной программе использованы
такие функции: expand - раскрытие скобок; collect - приведение подобных членов по
степеням указанной переменной; simplify - упрощение выражения; subs - подстановка.
Например, первый оператор выполняет подстановку
, (4.9)
а второй -
, .
Функция eval вычисляет значение символьного выражения при заданных значениях
его переменных, а функция sym2poly преобразует символьный полином в вектор его
коэффициентов.
В результате выполнения этой программы на экран будет выведена следующая
информация:
P =
[ w*(a+i*b), w*(a-i*b)]
G =
z^2-2*z*exp(T*w*a)*cos(T*w*b)+exp(2*T*w*a)
G_bin =
z^2-2*z*exp(-T*w)+exp(-2*T*w)
T =
0.0050
w =
200
G_bin_c =
1.0000 -0.7358 0.1353
Таким образом, мы получили, что при комплексно-сопряженных корнях (4.8)
аналогового полинома 2-го порядка дискретный полином, полученный методом
соответствия корней (см. выражения (4.6), (4.7)), имеет вид:
; (4.10)
Дискретный полином, соответствующий аналоговому прототипу с биномиальными
коэффициентами, имеющему одинаковые действительные корни , выглядит так:
(4.11)
и, наконец, этот же полином при , принимает вид
.
Подставляя в (10) , получим выражение для дискретной аппроксимации
полинома Баттерворта 2-го порядка.
Конечно же, все преобразования для полиномов второго порядка достаточно
легко выполнить и вручную. Однако с повышением степени полинома значительно
возрастает как сложност ь преобразований, так и количество возможных сочетаний
корней аналогового прототипа.
Приведем выражения для дискретных многочленов третьего-пятого порядков,
полученные при помощи программ с использованием средств символьной математики,
аналогичных приведенной выше.
1)
а),:
;
б), - полином Баттерворта:
;
в) - полином с биномиальными коэффициентами:
.
2)
а) , :
,
где
, , ;
б) , - кратные комплексно-сопряженные корни (при
- полином, сконструированный методом двойных пропорций):
,
где
,;
в) , - полином с биномиальными коэффициентами:
;
3)
а) ,, :
,
где
,,;;
б),,:
,
где
,,;
в), - полином с биномиальными коэффициентами:
.
Анализируя полученные результаты, можно сделать некоторые обобщения.
Прежде всего, используя метод математической индукции, можно записать
формулу дискретной аппроксимации полинома n-ой степени с биномиальными
коэффициентами:
, (4.12)
где
(4.13)
- биномиальные коэффициенты нормированного по СГК аналогового прототипа.
Коэффициенты бинома Ньютона (4.13) k-го порядка можно получить сложением
двух расположенных рядом биномиальных коэффициента полинома (k-1)-го порядка:
, (4.14)
для чего используют схему, называемую треугольником Паскаля:
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
1 6 15 20 15 6 1
1 7 21 35 35 21 7 1
1 8 28 56 70 56 28 8 1
...
Для дополнительной проверки на наличие ошибок, допускаемых по
невнимательности, стоит отметить, что при любом составе корней аналогового
полинома для коэффициентов разложения (4.7) справедливы следующие
утверждения:
1) сумма модулей коэффициентов при выражениях, относящихся к , равна
соответствующему биномиальному коэффициенту ;
2) сумма модулей всех коэффициентов при выражениях этого разложения равна 2n ;
3) суммы модулей всех коэффициентов при выражениях, относящихся к четным и
нечетным степенях z, равны друг другу и составляют 2n-1 ;
4) свободный член дискретной передаточной функции (4.7) и ее коэффициент при zn-1
можно определить по формулам:
, (4.15)
, (4.16)
где
;;;
, - действительные и мнимые части комплексно-сопряженных пар полюсов аналогового
прототипа;
- его действительные полюсы.
Формулы для расчета других коэффициентов дискретного полинома будут
отличаться друг от друга в зависимости от порядка полинома и состава корней
аналогового прототипа.
Из общих соотношений приведем еще одну формулу, справедливую для дискретных
полиномов четного порядка, аналоговые прототипы которых имеют только кратные
комплексно-сопряженные корни (,):
, (4.17)
. (4.18)
Приведенные положения и формулы можно использовать для контроля результатов
выполнения дискретной аппроксимации в аналитической форме аналоговых полиномов
более высоких порядков.
Полученные по предлагаемой методике коэффициенты дискретного полинома или
выражения для них могут быть в дальнейшем использованы для синтеза регуляторов и
наблюдателей состояния полного порядка. Методика автоматизации этого процесса в
среде пакета MatLab совпадает с изложенной ранее методикой синтеза аналоговых
регуляторов и наблюдателей состояния [14].
ВЫВОДЫ
1. Задача автоматизации формирования дискретного полинома по корням его
аналогового прототипа достаточно просто решается в среде пакета MATLAB как в
численном виде, так и в аналитической форме.
2. Дискретную аппроксимацию аналогового полинома с биномиальными коэффициентами
(или, что одно и тоже, с кратными действительными корнями) можно выполнять,
используя формулы (4.12)-(4.14).
3. Для контроля правильности выполненной дискретной аппроксимации можно
использовать специфические свойства дискретных полиномов, сформулированные
в данной статье, и формулы (4.16)-(4.18).
ЛИТЕРАТУРА
1. Куо Б. Теория и проектирование цифровых систем управления. -
М.: Машиностроение, 1986. 448 с.
2. Голубь А.П., Кузнецов Б.И., Опрышко И.А., Соляник В.П. Системы управления
электроприводами: Учеб. пособие. Киев: УМК ВО. 1992. 376с.
3. Kailath, T. Linear Systems. Prentice-Hall, 1980.
4. Kautsky, Nichols, Van Dooren. Robust Pole Assignment in Linear State
Feedback // Intl. J. Control, 41(1985)5. Р. 1129-1155
5. Медведев В.С., Потемкин В.Г. Control System Toolbox. МАТЛАБ 5 для студентов.
М.: ДИАЛОГ-МИФИ, 1989. 287 с.
6. Butterworth S. On the theory of filter amplifiers // Wireless Engineer,
London. - 1930. Vol. 7.
7. Айзинов М. М. Анализ и синтез линейных радиотехнических цепей в переходном
режиме. - Л.: Энергия, 1968. - 376с.
8. Кузовков Н.Т. Модальное управление и наблюдающие устройства. -
М.: Машиностроение, 1976. - 184 с.
9. Лазарев Ю.Ф. MatLab 5.x. Киев: BHV, 2000. 384 с.
10. Дьяконов В., Круглов В. Математические пакеты расширения MATLAB: Специальный
справочник. СПб.: Питер, 2001. 480с.
11. Толочко О.И., Коцегуб П.Х., Розкаряка П.И. Конструирование передаточных
функций линейных САУ по заданному перерегулированию // Вісник Національного
технічного університету "Харківський політехнічний інститут". Збірка наукових
праць. Тематичний випуск 10. - Харків: НТУ ХПІ, 2001. - С. 95-98.
12. D. Graham. R. Lathrop. The Synthesis of Optimum Transients Response. Criteria
and Standard Form. // AJEE, tech. - 1953. - P. 53-249.
13. Толочко О. И., Коцегуб П. Х., Губарь Ю. В., Федоряк Р. В. Конструирование
передаточных функций линейных САУ из условий модульного оптимума // Збiрник
наукових праць ДонДТУ, вип. 17: Донецьк: ДонДТУ. - 2000. - С. 24-30.
14. Толочко О.И., Федоряк Р.В. Автоматизация синтеза регуляторов и наблюдателей
состояния в среде пакета MatLab // Труды всероссийской научной конференции
" Проектирование научных и инженерных приложений в среде MatLab "
(ISBN 5-201-14940-5) - М.: ИПУ РАН, 2002. с. 482-496.
а) б)
Рис. 1.1. Структурные схемы САУ с регулятором состояния полного порядка
а) б)
Рис. 1.2. Структурные схемы наблюдателя состояния полного порядка
Синтез РСПП и НСПП состоит в выборе коэффициентов корректирующих обратных
связей K и L соответственно из условия обеспечения желаемого распределения
полюсов или, что одно и то же, желаемых коэффициентов характеристического
полинома замкнутой САУ с РСПП или замкнутого НСПП.
Таким образом, для определения искомых коэффициентов необходимо решить одно
из уравнений
(1.1)
или
, (1.2)
где Aз - матрица состояния замкнутой системы, которая определяется для САУ
с РСПП как
для НСПП как
- вектор желаемых полюсов замкнутой системы;
- желаемый характеристический полином замкнутой системы.
Существуют также аналитические выражения для вычисления искомых матриц
коэффициентов РСПП и НСПП. Так, в соответствии с опубликованной в 1980г. формулой
Аккермана [3] коэффициенты РСПП можно определить из последней строки матрицы
(1.3)
где
- матрица управляемости.
Ещё один способ основывается на использовании прямого и обратного
преобразования математического описания разомкнутой и замкнутой систем в
пространстве состояний к канонической форме Фробениуса. Согласно алгоритму
синтеза, приведенному в [1, 2], уравнения
и
имеют аналитические решения:
(1.4)
(1.5)
где
- коэффициенты характеристического полинома разомкнутой системы (или разомкнутого
наблюдателя);
- коэффициенты желаемого характеристического полинома;
- матрица наблюдаемости;
, 
.
Другой способ поиска вектора искомых коэффициентов требует представления
математического описания системы управления в канонической форме фазовой
переменной, или перехода к описанию в канонической форме, если исходные уравнения
состояния записаны в естественном базисе [1]. Пара матриц А и B считается
записанной в каноническом базисе, если эти матрицы имеют вид:
, 
Следует отметить, что пара матриц А и А* в естественном и каноническом базисах
имеют одни и те же собственные значения, кроме того, существует матрица
преобразования (перехода) к канонической форме Р, такая, что выполняются условия:
, 
Неособая матрица Р при этом задаётся в виде:
где
Представив описание объекта в канонической форме фазовой переменной,
вычисляют вектор коэффициентов регулятора состояния в каноническом базисе как
(1.6)
где 
- коэффициенты характеристических полиномов соответственно объекта
(разомкнутой системы) и желаемого (замкнутой системы).
Имея матрицу перехода к каноническому базису Р, можно получить действительные
коэффициенты обратных связей
(1.7)
Методики, описываемые формулами (1.6, 1.7) и (1.4, 1.5), удобны тем, что
в выражения входят разности коэффициентов при соответствующих степенях
характеристических полиномов разомкнутой системы и желаемого стандартного
полинома. А все остальные задачи - формирование матриц и матричные операции,
предусмотренные тем или иным методом синтеза - хорошо поддаются алгоритмизации,
и не представляют затруднений.
Формулы (1.1) и (1.2), (1.3), (1.4) и (1.5), (1.6) и (1.7) могут с успехом
лечь в основу как численных, так и аналитических методов автоматизации синтеза
регуляторов и наблюдателей состояния.
Используемые в аналитических методах чинтеза РСПП и НСПП матрицы
преобразования при высоком порядке дифференциальных уравнений, описывающих объект
регулирования, могут оказаться плохо предопределенными и давать большую
погрешность. В связи с этим заслуживают внимания также и численные алгоритмы
синтеза регуляторов и наблюдателей состояния полного порядка, наиболее популярный
из которых предложен в [4].
Следует отметить, что, если известен какой-либо алгоритм поиска коэффициентов
РСПП, то подставляя в него вместо матриц A и B матрицы AТ и CТ соответственно,
можно определить транспонированный вектор коэффициентов НСПП LТ. Эта особенность
обусловлена дуальной связью между управляемостью и наблюдаемостью.
Итак, при описании алгоритмов синтеза РСПП и НСПП мы будем рассматривать
четыре основных способа:
1. Алгоритм, основанный на составлении и решении системы n линейных
алгебраических уравнений с n неизвестными относительно искомых коэффициентов
(в матричной трактовке этот алгоритм описывается формулой (1.2)).
2. Алгоритмы расчета, использующие стандартные (т.е. входящие в стандартную
поставку системы MatLab) функции. Это прежде всего, способ, основанный на
использовании формулы Аккермана (1.3), и численный метод, предложенный в [4].
3. Алгоритм, описанный в [1,2] и основанный на использовании формул (1.4) и (1.5).
4. Алгоритм вычислений с использованием перехода к описанию объекта в
канонической форме фазовой переменной (формулы (1.6),(1.7)).
а)
б)
Рис. 2.1 Структурные схемы объекта "Тиристорный преобразователь - двигатель
постоянного тока".
При реальном замыкании обратных связей после их расчета, следует определить
очередность следования переменных состояния при помощи команды
[[sizes,x0,xstr]=имя_файла_модели
Интересующая нас информация будет содержаться в переменной xstr в виде путей
к переменным состояния, например,
xstr =
'TPD/Integrator1'
'TPD/Integrator2'
'TPD/Integrator'
некоторого аналогового полинома n-ой степени
, нормированного по его среднегеометрическому корню (СГК)
:
, (4.1)
, (4.2)
, (4.3)
. (4.4)
Один из наиболее удобных способов дискретизации полинома (4.1) состоит из
следующих этапов:
1) определение корней 
полинома (4.3) и корней 
полинома (1):
; (4.5)
2) определение корней дискретного полинома:
, (4.6)
где Т - период дискретности;
3) формирование дискретного полинома
. (4.7)
В пакете MatLab информация о полиномах хранится в виде массива
коэффициентов, упорядоченных по убыванию степени независимой переменной. Например,
вектор
содержит информацию о полиноме (4.3).
Корни полинома по его коэффициентам определяет функция roots, например,
P=roots(G), а обратную операцию выполняет функция poly, например, G=real(poly(P)).
Следует отметить, что в справочной литературе могут быть даны как
коэффициенты, так и корни стандартных полиномов, однако чаще приводятся значения
коэффициентов . Поэтому сведем в табл. 1 значения нормализованных корней
стандартных полиномов, наиболее часто используемых при синтезе систем управления
электромеханическими объектами: 1 - полином Баттерворта [6-8], 2 - полином
с 5%-ным перерегулированием [11], 3 - полином Грехема-Летропа [12, 8], 4 - полином
Бесселя [7], 5 - полином, сформированный методом двойных пропорций [13],
6 - полином с кратными комплексно-сопряженными корнями и коэффициентом
демпфирования 0,75. В таблице отсутствуют корни полиномов с биномиальными
коэффициентами, поскольку все они равны -1.
определяются следующими
выражениями:
,
.
Несколько сложнее поставленная задача решается в символьном виде. В этом
случае нам необходимо найти не численные значения коэффициентов дискретного
полинома, а аналитические зависимости этих коэффициентов от параметров корней его
аналогового прототипа. Решение задачи усложняется тем, что вид указанных
аналитических выражений зависит от того, сколько корней образуют комплексно-
сопряженные пары и сколько из них являются действительными, а также от
соотношения амплитуд и фаз, действительных и мнимых составляющих этих корней.
Вариант программы, решающей поставленную задачу для полиномов 2-го порядка,
может иметь такой вид:
syms z p a b T
P(1)=a+i*b; P(2)=w*a-i*b
Z=exp(T*P);
G=prod(z-Z);
G=expand(G);
G=subs(G,exp(i*T*b),cos(T*b)+i*sin(T*b));
G=simplify(G);
G=collect(G,z)
G_bin=subs(G,[a b],[-w 0])
T=0.005, w=1/T
G_bin_c=sym2poly(eval(G_bin))
Описатель syms объявляет следующие за ним переменные символьными.
Для работы с комплексными числами в пакете MatLab предопределены
переменные i и j, значением которых является 
. Таким образом, в программе
первоначально заданы два комплексно-сопряженных корня
. (4.8)
Функция prod вычисляет произведение всех элементов вектора-аргумента. Из
функций преобразования символьных выражений в данной программе использованы
такие функции: expand - раскрытие скобок; collect - приведение подобных членов по
степеням указанной переменной; simplify - упрощение выражения; subs - подстановка.
Например, первый оператор выполняет подстановку
, (4.9)
а второй -
,
.
Функция eval вычисляет значение символьного выражения при заданных значениях
его переменных, а функция sym2poly преобразует символьный полином в вектор его
коэффициентов.
В результате выполнения этой программы на экран будет выведена следующая
информация:
P =
[ w*(a+i*b), w*(a-i*b)]
G =
z^2-2*z*exp(T*w*a)*cos(T*w*b)+exp(2*T*w*a)
G_bin =
z^2-2*z*exp(-T*w)+exp(-2*T*w)
T =
0.0050
w =
200
G_bin_c =
1.0000 -0.7358 0.1353
Таким образом, мы получили, что при комплексно-сопряженных корнях (4.8)
аналогового полинома 2-го порядка дискретный полином, полученный методом
соответствия корней (см. выражения (4.6), (4.7)), имеет вид:
; (4.10)
Дискретный полином, соответствующий аналоговому прототипу с биномиальными
коэффициентами, имеющему одинаковые действительные корни 
(4.11)
и, наконец, этот же полином при 
,
принимает вид
.
Подставляя в (10) 
, получим выражение для дискретной аппроксимации
полинома Баттерворта 2-го порядка.
Конечно же, все преобразования для полиномов второго порядка достаточно
легко выполнить и вручную. Однако с повышением степени полинома значительно
возрастает как сложност ь преобразований, так и количество возможных сочетаний
корней аналогового прототипа.
Приведем выражения для дискретных многочленов третьего-пятого порядков,
полученные при помощи программ с использованием средств символьной математики,
аналогичных приведенной выше.
1)
а)
,
:
;
б)
,
- полином Баттерворта:
;
в)
- полином с биномиальными коэффициентами:
.
2)
а)
,
:
,
где
,
,
;
б)
,
- кратные комплексно-сопряженные корни (при
- полином, сконструированный методом двойных пропорций):
,
где
,
;
в) 
,
- полином с биномиальными коэффициентами:
;
3)
а) 
,
,
:
,
где
,
,
;
;
б)
,
,
:
,
где
,
,
;
в)
,
- полином с биномиальными коэффициентами:
.
Анализируя полученные результаты, можно сделать некоторые обобщения.
Прежде всего, используя метод математической индукции, можно записать
формулу дискретной аппроксимации полинома n-ой степени с биномиальными
коэффициентами:
, (4.12)
где
(4.13)
- биномиальные коэффициенты нормированного по СГК аналогового прототипа.
Коэффициенты бинома Ньютона (4.13) k-го порядка можно получить сложением
двух расположенных рядом биномиальных коэффициента полинома (k-1)-го порядка:
, (4.14)
для чего используют схему, называемую треугольником Паскаля:
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
1 6 15 20 15 6 1
1 7 21 35 35 21 7 1
1 8 28 56 70 56 28 8 1
...
Для дополнительной проверки на наличие ошибок, допускаемых по
невнимательности, стоит отметить, что при любом составе корней аналогового
полинома для коэффициентов разложения 
(4.7) справедливы следующие
утверждения:
1) сумма модулей коэффициентов при выражениях, относящихся к 
, равна
соответствующему биномиальному коэффициенту 
;
2) сумма модулей всех коэффициентов при выражениях этого разложения равна 2n ;
3) суммы модулей всех коэффициентов при выражениях, относящихся к четным и
нечетным степенях z, равны друг другу и составляют 2n-1 ;
4) свободный член дискретной передаточной функции (4.7) и ее коэффициент при zn-1
можно определить по формулам:
, (4.15)
, (4.16)
где
;
;
;
,
- действительные и мнимые части комплексно-сопряженных пар полюсов аналогового
прототипа;
- его действительные полюсы.
Формулы для расчета других коэффициентов дискретного полинома будут
отличаться друг от друга в зависимости от порядка полинома и состава корней
аналогового прототипа.
Из общих соотношений приведем еще одну формулу, справедливую для дискретных
полиномов четного порядка, аналоговые прототипы которых имеют только кратные
комплексно-сопряженные корни (
,
):
, (4.17)
. (4.18)
Приведенные положения и формулы можно использовать для контроля результатов
выполнения дискретной аппроксимации в аналитической форме аналоговых полиномов
более высоких порядков.
Полученные по предлагаемой методике коэффициенты дискретного полинома или
выражения для них могут быть в дальнейшем использованы для синтеза регуляторов и
наблюдателей состояния полного порядка. Методика автоматизации этого процесса в
среде пакета MatLab совпадает с изложенной ранее методикой синтеза аналоговых
регуляторов и наблюдателей состояния [14].