Матрица Мироновского. Грамиан M=S(.)-1S(.) на основе двух матриц Сильвестра, отличающихся аргументами. Позволяет вычислять ганкелевы числа динамической системы непосредственно по коэффициентам полиномов числителя и знаменателя передаточной функции.

Рассмотрим реакцию динамической системы с передаточной функцией

Q(p)=b(p)/a(p)

на входной сигнал u(p)=c(-p)/a(-p), который зеркально-симметричен собственному движению системы и отличается масштабом y(p)=sc(p)/a(p). Полином a(p) - гурвицев. Эта реакция имеет вид

Q(p)u(p)=b(p)c(-p)/a(p)a(-p)=sc(p)/a(p)+d(p)/a(-p),

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

Задача состоит в том, чтобы найти неизвестные заранее ганкелевы числа (коэффициенты передачи s) и ганкелевы функции (соответствующие полиномы с(p)). Для этого рассмотрим полиномиальное уравнение b(p)c(-p)=sa(-p)c(p)+a(p)d(p), следующие из формулы выше. Перестановкой членов получим

s(a(p) a(-p))(d(p)/s c(p))'=(0 b(p))(d(-p)/s c(-p))',

что приводит к матричному уравнению sS(a(p),a(-p))x=S(0,b(p))Zx, где S - матрица Сильвестра, Z=diag(1 -1.. 1 -1..) - сигнатурная матрица, согласующая смену знака аргументов x, которые содержит искомые коэффициенты d(p)/s и c(p). Усечением вектора x до колонки коэффициентов с(p), задача сводится к стандартной полной проблеме собственных значений

Mx=sx,

где M=[S(a(p),a(-p))-1]S(b(p))Z - кроcс-грамиан (матрица Мироновского), [.] - проектор, оставляющий лишь нижние блоки, здесь Z - усеченная сигнатурная матрица, чередующая знаки столбцов.

Можно составить уравнения иначе

(a(p) b(p))(-d(p) c(-p))'=s(0 a(-p))(-d(-p) c(p))',

что приводит к матричному уравнению S(a(p),b(p))x=sS(0,a(-p))Zx, где x содержит коэффициенты -d(p) и c(-p). Инверсией S(a(p),b(p)) и усечением вектора x до колонки коэффициентов с(-p), задача сводится к стандартной полной проблеме собственных значений с уравнением вида

Mx=kx,

где кросс-грамиан M=[S(a(p),b(p))-1]S(a(-p))Z. Искомые значения s инверсны собственным значениям k, а x получаем как собственные векторы. Эти процедуры вошли в системный тулбокс control, обращения к которому позволяют вычислить ганкелевы числа и функции.


%toolbox control, 
b=[1 2], a=[1 2 1], t=time(10),

% НАХОЖДЕНИЕ СОБСТВЕННЫХ ВЕКТОРОВ,
af=flip(a), bf=flip(b), 
M=mironovsky(af bf), 
[c s]=eig(M),  c=flips(c), 
diag(s), s=?,   

% ПРОВЕРКА (ВТОРОЙ СПОСОБ),
h=heig(b a), c2=h[0], s2=h[1], s2=?, 
 
% РАСЧЕТ ГАНКЕЛЕВЫХ ФУНКЦИЙ,
f1=impulse(c[0] a t), 
f2=impulse(c[1] a t), plot(t [f1 f2]),

Обратная задача. Что меняется в задаче, если, наоборот, a(p),c(p) известны (пусть и s известно)? Возникает линейное уравнение a(p)(-d(p))+c(-p)b(p)=sa(-p)c(p) относительно -d(p),b(p), в котором нас интересует только b(p), порождающее матричное

b=sMc,

где M=[S(a(p),c(-p))-1]S(a(-p)) - отличается от кросc-грамиана предыдущей задачи аргументами и сигнатурой. Результант det(S(a(p),c(-p))) позволяет оценивать разрешимость. Возьмем в качестве c одно из решений выше, получим b.


%toolbox control, 
a=[1 2 1], s=1.059, c=[0.5257 0.85],  

% НАХОЖДЕНИЕ ПЕРЕДАТОЧНОЙ ФУНКЦИИ ПО ГАНКЕЛЕВОЙ ФУНКЦИИ,
af=flip(a), cf=flip(c), cz=signature(cf), 
M=crossmatrix(af cz), M=signature(M), 
b={M*cf}, b={s*b}, b=flip(b), b=?,

% ВТОРИЧНОЕ ВЫЧИСЛЕНИЕ ГАНКЕЛЕВОЙ ФУНКЦИИ,
t=time(10), b=[1 2], 
f=impulse(c a t), 
% ПРОВЕРКА ГАНКЕЛЕВЫМ ЭКСПЕРИМЕНТОМ,
y=hsim(b a f t),  plot(t [f y]),

Примечание. Для фазовращателя c Q(p)=a(-p)/a(p)=(-p+1)(-p+2)/(p+1)(p+2) ганкелева функция набирается из 'кусочков': y(p)=(-p+1)/(p+2). Входной сигнал u(p)=(p+1)/(-p+2) воспроизводит эту функцию. Более того, это касается класса передаточных функций, отличающихся слагаемым (совпадение ганкелевых функций у разных подсистем, к тому же, разных порядков): Q(p)+k/(p+1).



Rambler's Top100