Бороденко В.А. Исследование систем управления в среде MATLAB
Подождите немного. Документ загружается.
171
Псевдоканоническая управляемая форма с сопровождающей
(companion) матрицей А, содержащей характеристический многочлен
в правом столбце (авторы MATLAB советуют избегать этой формы
ввиду численной неустойчивости) может быть получена обращением
>> sys=canon(w,'companion')
a =
x1 x2 x3
x1 0 0 -15
x2 1 0 -11
x3 0 1 -5
b =
u1
x1 1
x2 0
x3 0
c =
x1 x2 x3
y1 1 -2 1
d =
u1
y1 0
MATLAB осуществляет преобразование описания системы к ле-
стничной форме с разделением управляемых и неуправляемых (соот-
ветственно наблюдаемых и ненаблюдаемых) мод командами
[Ap,Bp,Cp,P,K]=ctrbf(A,B,C) и [Ap,Bp,Cp,P,K]=obsvf(A,B,C). Исполь-
зуя неособую матрицу преобразования подобия P так, что A
p
= PAP
T
,
B
p
= PB, C
p
= CP
T
, система приводится к форме
cncp
c
p
c
nc
p
CCC
B
B
AA
A
A
,
0
,
0
21
,
где пара (A
c
, B
c
) управляема, причем С
с
(s1-A
c
)
-1
B
c
= C(s1-A)
-1
B, либо к
форме
op
o
no
p
o
no
p
CC
B
B
B
A
AA
A 0 , ,
0
12
где пара (A
o
, B
o
) наблюдаема, причем С
o
(s1-A
o
)
-1
B
o
= C(s1-A)
-1
B. Эле-
менты вверху и слева с индексом nc описывают неуправляемую, а с
индексом no – ненаблюдаемую часть системы. Левый аргумент К
представляет собой вектор длины n и содержит число управляемых
(наблюдаемых) состояний, определяемых на каждом шаге вычисления
матрицы P. Предусмотрено задание погрешности преобразования tol в
форме ctrbf(A,B,C,tol) и obsvf(A,B,C,tol). В приводимом примере век-
172
тор k содержит три управляемых моды и лишь одну наблюдаемую -3.
>> a=[0 1 0;0 0 1;-6 -11 -6]; b=[0 0 1]'; c=[2 3 1]; d=0;
>> sys=ss(a,b,c,d);
>> w=tf(sys)
Transfer function: % ПФ исходной системы
s^2 + 3 s + 2
----------------------
s^3 + 6 s^2 + 11 s + 6
>> rank(ctrb(sys)) % система полностью управляема
ans =
3
>> [ap,bp,cp,p,k]=ctrbf(a,b,c) % разделяем моды
ap =
0 1 0
0 0 1
-6 -11 -6
bp =
0
0
1
cp =
2 3 1
p = % система не изменилась
1 0 0
0 1 0
0 0 1
k =
1 1 1 % три управляемых моды
>> rank(obsv(sys)) % система частично наблюдаема
ans =
1
>> [ap,bp,cp,p,k]=obsvf(a,b,c)
ap = % плохо вычислены ненаблюдаемые моды -2 и -1
-1.8833 1.9244 -11.4475
0.0536 -1.1167 6.8158
-0.0000 -0.0000 -3.0000 % управляемая мода -3
bp =
0.8465
-0.4604
0.2673
cp =
0.0000 0.0000 3.7417
p =
0.2673 -0.4604 0.8465
0.8018 -0.3811 -0.4604
0.5345 0.8018 0.2673
173
k =
1 0 0 % одна управляемая мода
>> syspo=ss(ap,bp,cp,0);
>> tf(syspo) % ПФ преобразованной системы не изменилась
Transfer function:
s^2 + 3 s + 2
----------------------
s^3 + 6 s^2 + 11 s + 6
Все формы модели в пространстве состояния, создаваемые
MATLAB, не совпадают с общепринятыми в СНГ каноническими
управляемым и наблюдаемым представлениями, основанными на ме-
тоде фазовых переменных. Этот пробел призваны компенсировать
специально разработанные и приводимые далее функции пользовате-
ля tf2cs() и tf2os(). Первая по числителю и знаменателю передаточной
функции выводит каноническую управляемую форму, причем предла-
гаются два варианта этой функции: с вычислением матриц делением
полинома числителя ПФ на полином знаменателя
function [a,b,c,d]=tf2cs(num,den)
a=flipud(fliplr(compan(den)));
b=flipud(eye(length(den)-1,1));
[d,c]=deconv(num/den(1),den/den(1));
c=fliplr(c); c=[c zeros(1,length(b))]; c=c(1:length(b));
и с вычислением элементов матрицы с по правилу c
i
= b
n+1-i
– a
n+1-i
∙d
function [a,b,c,d]=tf2cs(num,den)
a=flipud(fliplr(compan(den)));
b=fliplr(eye(1,length(den)-1))';
num=[zeros(1,length(den)-length(num)) num/den(1)];
d=num(1); c=num-den/den(1)*d; c=fliplr(c(2:end));
Результат работы обоих вариантов функции совпадает до нуле-
вого порядка числителя и знаменателя передаточной функции
>> num=[1 3 2]; den=[1 5 11 15];
>> [a, b, c, d]=tf2cs(num,den)
a =
0 1 0
0 0 1
-15 -11 -5
b =
0
0
174
1
c =
2 3 1
d =
0
Вторая функция tf2os()
function [a,b,c,d]=tf2оs(num,den)
a=flipud(fliplr(compan(den)));
ld=length(den); c=eye(1,ld-1);
num=[zeros(1,ld-length(num)) num/den(1)]; den=den/den(1);
b(1)=num(1); n=1;
while n<ld
n=n+1;
b(n)=num(n)-sum(b.*fliplr(den(2:n)));
end
d=b(1); b=b(2:ld)';
по числителю и знаменателю передаточной функции возвращает тра-
диционное описание системы в канонической наблюдаемой форме
>> [a, b, c, d]=tf2оs(num,den)
a =
0 1 0
0 0 1
-15 -11 -5
b =
1
-2
1
c =
1 0 0
d =
0
MATLAB позволяет упростить модель системы регулирования
(снизить ее порядок), используя функции:
modred(sys,elim) – снижает порядок модели, описанной в про-
странстве состояний, исключая состояния, обозначенные вектором
индексов elim. Последний можно сформировать с применением функ-
ции balreal(), выявляющей переменные с пренебрежимым вкладом в
реакцию системы.
>> [sys,g] = balreal(sys) % вычисление сбалансированной реализации
>> elim = (g<1e-8) % индикация пренебрежимых состояний
>> rsys = modred(sys,elim) % удаление пренебрежимых состояний
175
minreal() – строит минимальную реализацию системы путем ис-
ключения неуправляемых и ненаблюдаемых переменных состояния
(пар нуль-полюс с одинаковыми значениями);
sminreal() – удаляет состояния, не имеющие связи ни с одним
входом или выходом;
balred() – возвращает аппроксимирующую модель пониженного
порядка, задаваемые свойства описывает balredOptions.
Преобразовывать модель системы из одного представления в
другое приходится в разнообразных целях, в частности, при проекти-
ровании модальных регуляторов необходимо перейти к канонической
форме управляемости и затем вернуться к исходному виду. Если обо-
значить новый вектор состояния как h(t), то матрицы коэффициентов
новой системы вычисляются по формулам A
h
=PAP
-1
, B
h
=PB, C
h
=CP
-1
,
k
h
= kP
-1
. Задаваясь произвольной неособой матрицей Р необходимого
размера, можно получить бесконечное множество описаний одной и
той же системы в пространстве состояний. Однако при любых преоб-
разованиях исходная и преобразованная система должны иметь оди-
наковые собственные значения (корни характеристического полино-
ма), а преобразование подобия не должно менять передаточную
функцию системы.
Рассмотрим преобразование произвольной системы к канониче-
ской управляемой форме, канонической наблюдаемой форме, запи-
санных методом фазовых переменных, и представлению с диагональ-
ной матрицей А, составленной из собственных значений системы, ис-
следуя вопрос получения матрицы преобразования Р.
Создадим сначала модель системы в пространстве состояний с
произвольным видом матриц A, B, C функцией ss().
>> w=tf([1 2 3],[1 6 11 6]);
>> sys=ss(w);
>> [A,b,c,d]=ssdata(sys)
A =
-6.0000 -2.7500 -1.5000
4.0000 0 0
0 1.0000 0
b =
2
0
0
c =
0.5000 0.2500 0.3750
d =
0
176
Для описания системы диагонализированной матрицей А, мат-
рицу Т собственных векторов системы, обратную к матрице преобра-
зования Р, и преобразованную матрицу Аd с комплексными значе-
ниями собственных чисел (корней характеристического полинома),
можно получить стандартной функцией eig().
>> [T,Ad]=eig(A)
T =
-0.5797 0.4082 0.1741
0.7730 -0.8165 -0.6963
-0.2577 0.4082 0.6963
Ad =
-3.0000 0 0
0 -2.0000 0
0 0 -1.0000
Используя функцию ss2ss(), создадим преобразованное к диаго-
нальному виду описание системы (по диагонали матрицы А размеща-
ются полюсы системы, в том числе комплексные)
>> P=inv(T); sysd=ss2ss(sys,P)
a =
x1 x2 x3
x1 -3 -8.438e-015 -7.994e-015
x2 1.155e-014 -2 -1.155e-014
x3 -3.109e-015 3.553e-015 -1
b =
u1
x1 -15.52
x2 -19.6
x3 5.745
c =
x1 x2 x3
y1 -0.1932 0.1531 0.1741
d =
u1
y1 0
проверим его правильность, восстановив передаточную функцию
>> tf(sysd)
Transfer function:
s^2 + 2 s + 3
----------------------
s^3 + 6 s^2 + 11 s + 6
и очистим матрицу А от величин малого порядка
177
>> sysd.a=sysd.a.*(abs(sysd.a)>1.e-10)
Рассмотрим приведение к канонической управляемой форме. В
этом случае вид матрицы А (сопровождающей или матрицы Фробе-
ниуса) полностью определяется характеристическим полиномом, вид
матрицы В также формализован. Таким образом, формируем матрицы
Ас и bс системы в канонической управляемой форме и вычисляем
матрицы управляемости обеих систем, тогда матрица преобразования
Р равна отношению матрицы управляемости новой системы к матрице
управляемости исходной.
>> Ac=fliplr(flipud(compan(poly(A))))
Ac =
0 1.0000 0
0 0 1.0000
-6.0000 -11.0000 -6.0000
>> bc=[0; 0; 1]; q=ctrb(A,b); qc=ctrb(Ac,bc); p=qc/q
p =
0 0 0.1250
0 0.1250 0
0.5000 0 0
>> sysc=ss2ss(sys,p)
a =
x1 x2 x3
x1 0 1 0
x2 0 0 1
x3 -6 -11 -6
b =
u1
x1 0
x2 0
x3 1
c =
x1 x2 x3
y1 3 2 1
d =
u1
y1 0
Проверим результат самостоятельными преобразованиями
>> Ac=p*A/p
Ac =
0 1 0
0 0 1
-6 -11 -6
>> cc=c/p
cc =
3 2 1
178
Обязательным условием преобразования является невырож-
денность матрицы управляемости исходной системы.
Процедура перехода к канонической наблюдаемой форме анало-
гична и отличается лишь тем, что используются матрицы наблюдае-
мости, причем матрица преобразования базиса вычисляется по отно-
шению матрицы наблюдаемости исходной системы к матрице наблю-
даемости новой. Очевидным условием возможности преобразования
является полная наблюдаемость новой системы.
>> Ao=fliplr(flipud(compan(poly(A))))
Ao =
0 1.0000 0
0 0 1.0000
-6.0000 -11.0000 -6.0000
>> co=[1 0 0]; n=obsv(A,c); no=obsv(Ao,co); P=n/no
P =
0.5000 0.2500 0.3750
-2.0000 -1.0000 -0.7500
8.0000 4.7500 3.0000
>> syso=ss2ss(sys,P)
a =
x1 x2 x3
x1 0 1 0
x2 0 0 1
x3 -6 -11 -6
b =
u1
x1 1
x2 -4
x3 16
c =
x1 x2 x3
y1 1 0 0
d =
u1
y1 0
Произведем проверку самостоятельным вычислением матриц А,
b новой системы и ее передаточной функции – результаты совпадают.
>> Aо=P*A/P
Aо =
0 1 0
0 0 1
-6 -11 -6
>> bo=P*b
bo =
179
1
-4
16
>> tf(syso)
Transfer function:
s^2 + 2 s + 3
----------------------
s^3 + 6 s^2 + 11 s + 6
Обратный переход (возвращение к исходной системе), напри-
мер, после выбора параметров модального регулятора, во всех случаях
осуществляется применением матрицы Р в обратном порядке, т. е.
A = P
-1
A
h
P, B=P
-1
B
h
, C=C
h
P, k = k
h
P или обратной матрицы inv(P) в
функции ss2ss(). Рассмотрим обратный переход для системы, приве-
денной к канонической наблюдаемой форме
>> sys=ss2ss(syso,inv(P))
a =
x1 x2 x3
x1 -6 -2.75 -1.5
x2 4 -1.066e-014 -7.105e-015
x3 0 1 0
b =
u1
x1 2
x2 0
x3 0
c =
x1 x2 x3
y1 0.5 0.25 0.375
d =
u1
y1 0
Функция xperm(sys,P) переупорядочивает состояния SS модели
согласно вектору перестановок P с новыми позициями переменных.
>> [a,b,c,d]=linmod('standard3');
>> sys=ss(a,b,c,d)
a =
x1 x2 x3
x1 0 0 1
x2 -4 -2 -3
x3 0 1 0
b =
u1
x1 0
x2 1
180
x3 0
c =
x1 x2 x3
y1 1 0 0
d =
u1
y1 0
Continuous-time model.
>> P = [3 1 2];
>> sys = xperm(sys, P)
a =
x1 x2 x3
x1 0 0 1
x2 1 0 0
x3 -3 -4 -2
b =
u1
x1 0
x2 0
x3 1
c =
x1 x2 x3
y1 0 1 0
d =
u1
y1 0
Continuous-time model.
Она, например, позволяет отсортировать переменные состояния
по алфавиту, если им присвоены имена.
>> set(sys,'StateName',{'x3', 'x1', 'x2'});
>> [x, P] = sort(sys.StateName);
>> sys = xperm(sys, P)
В заключение напомним о том, что SS, TF и ZPK модели не яв-
ляются эквивалентными с точки зрения численных расчетов. Полюсы
передаточных функций, особенно высокого порядка, могут отличать-
ся от полюсов ZPK или SS моделей, из которых они были получены.
Преобразование в пространство состояний не является однозначным
для SISO системы и не гарантирует получения минимальной конфи-
гурации MIMO системы. При двойном преобразовании, например,
командой ss(tf(sys_ss)), может быть сформирована SS модель не толь-
ко с другими матрицами A, B, C, D, но даже с иным количеством пе-
ременных состояния. Поэтому желательно избегать излишних преоб-
разований форм и моделей, и нужно быть готовым к тому, что неожи-
данные для пользователя преобразования осуществит сама MATLAB.