Автоматика и телемеханика, № 2, 2019
Нелинейные системы
© 2019 г. А.Н. ЖИРАБОК, д-р техн. наук (zhirabok@mail.ru),
А.Е. ШУМСКИЙ, д-р техн. наук (a.e.shumsky@yandex.com)
(Дальневосточный федеральный университет, Владивосток,
Институт проблем морских технологий ДВО РАН, Владивосток)
НЕПАРАМЕТРИЧЕСКИЙ МЕТОД ДИАГНОСТИРОВАНИЯ
НЕЛИНЕЙНЫХ ДИНАМИЧЕСКИХ СИСТЕМ1
Рассматривается задача функционального диагностирования нелиней-
ных динамических систем непараметрическим методом, особенность ко-
торого состоит в том, что значения некоторых параметров диагностируе-
мой системы могут быть неизвестными. Предлагается подход, позволяю-
щий производить локализацию дефектов. Устанавливаются необходимые
и достаточные условия возможности преобразования системы к канони-
ческому виду, необходимому для применения непараметрического метода
при диагностировании заданной нелинейной системы.
Ключевые слова: нелинейные системы, диагностирование, поиск дефек-
тов, непараметрический метод.
DOI: 10.1134/S0005231019020028
1. Введение и постановка задачи
Настоящая статья является логическим продолжением работы [1], где рас-
сматривалась задача диагностирования технических систем, описываемых
динамическими линейными моделями, на основе непараметрического мето-
да, особенность которого состоит том, что знание значений постоянных па-
раметров системы не является обязательным для реализации процедуры ди-
агностирования, в линейном случае достаточно знать только размерности
векторов входа, состояния и выхода и места вхождения этих параметров в
описание системы [2-6].
В настоящей работе рассматриваются системы, описываемые нелинейны-
ми динамическими моделями. Напомним, что для применения непараметри-
ческого метода исходная система должна быть представлена в виде вход-
выходных соотношений, что всегда можно сделать для линейных моделей.
В отличие от этого далеко не всякая нелинейная система может быть пред-
ставлена в таком виде, поэтому основное внимание в работе уделяется разра-
ботке условий такого представления и процедур его получения.
Для формулирования этих условий и процедур используется специальный
математический аппарат - алгебра функций, разработанная в [7] для решения
1 Работа выполнена при финансовой поддержке Российского научного фонда (проект
№ 16-19-00046).
24
задач анализа нелинейных динамических систем и позволяющая рассмат-
ривать системы с негладкими нелинейностями, к числу которых относятся
такие типичные нелинейности, как сухое трение, люфт, гистерезис, насыще-
ние. Кроме этого, для решения задачи используется логико-динамический
подход [8], позволяющий при ограничении на класс используемых преобра-
зований исходной системы применять для решения задачи только методы
линейной алгебры.
Рассматривается класс систем, описываемых нелинейной динамической
моделью
(1.1)
x(t + 1) = f(x(t), u(t), γ), y(t) = h(x(t)).
Здесь x ∈ X ⊆ Rn, u ∈ Rm, y ∈ Rl - векторы состояния, управления и выхо-
да; γ = (γ1, . . . , γs)T - вектор параметров, i-й дефект ρi в системе проявляет-
ся в виде недопустимого отклонения i-го параметра γi от его номинального
значения, f и h - нелинейные функции, при этом функция f может быть
недифференцируемой.
Для решения задачи диагностирования непараметрическим методом тре-
буется произвести преобразование векторов состояния и выхода
(1.2)
x∗i(t) = ϕ(i)(x(t)), i = 1,... ,N, y
(t) = ψ(y(t))
с помощью некоторых функций ϕ(1), . . . , ϕ(N) и ψ такое, что в новых коорди-
натах преобразованная система описывается моделью без обратных связей:
x1(t + 1) =
(x2(t), . . . , x∗N (t), y(t), u(t)),
x∗i(t + 1) =
(x∗i+1(t), . . . , x∗N (t), y(t), u(t)), i = 2, . . . , N - 1,
(1.3)
x∗N(t + 1) =
(y(t), u(t)),
y(t) = h(x1(t))
для некоторых функций
,...,
, h; для упрощения будем опускать
вектор γ там, где это возможно. Отметим, что в нелинейном случае некоторые
переменные x∗i могут быть векторными.
В ряде случаев исходная модель может быть приведена к более простому
виду, нежели (1.3):
x1(t + 1) =
(x2(t), y(t), u(t)),
x∗i(t + 1) =
(x∗i+1(t), y(t), u(t)), i = 2, . . . , N - 1,
(1.4)
x∗N(t + 1) =
(y(t), u(t)),
y(t) = h(x1(t)).
Условия и способы приведения модели (1.1) к виду (1.3) или (1.4) рассмот-
рены в разделах 3 и 4. Остановимся вначале на использовании полученных
моделей для целей диагностирования.
25
2. Основные конструкции непараметрического метода
2.1. Преобразование построенной модели
Для применения непараметрического метода произведем в представлени-
ях (1.3) и (1.4) ряд временных сдвигов и подстановок, для конкретности рас-
смотрим только (1.4):
x1(t + 2) =
(x2(t + 1), y(t + 1), u(t + 1)) =
=
(
(x3(t), y(t), u(t)), y(t + 1), u(t + 1)),
x1(t+3) =
(
(
(x4(t), y(t), u(t)), y(t + 1), u(t + 1)), y(t + 2), u(t + 2)),
x1(t + N) = F(y(t),u(t),y(t + 1),u(t + 1),... ,y(t + N - 1),u(t + N - 1))
для некоторой функции F. В результате с учетом y = h(x1) получаем
(
(2.1) y(t + N) = F∗∗ y(t), u(t), y(t + 1), u(t + 1), . . .
)
...,y(t + N - 1),u(t + N - 1) ,
где F∗∗ = h(F).
Будем полагать, что функция F∗∗ представляет собой сумму нелинейных
членов вида
F∗∗ =
Γi(γ)Pi(y(t),u(t),y(t + 1),u(t + 1),... ,y(t + N - 1),u(t + N - 1)),
i=1
где Γi(γ) - алгебраическое выражение, представляющее собой функцию эле-
ментов вектора параметров γ, значения которых могут быть неизвестными,
функция Pi() может содержать только такие элементы вектора γ, значения
которых известны, i = 1, . . . , p, p - число слагаемых. Отметим, что искомое
представление всегда может быть получено для полиномиальных функций.
Если функция Pi() содержит параметр, значение которого может быть неиз-
вестно, то для возможности применения непараметрического метода ее сле-
дует представить степенным рядом, коэффициенты которого будут зависеть
от этого параметра.
Из сказанного следует, что равенство (2.1) может быть записано в виде
y(t + N) = F∗∗(y(t),u(t),... ,y(t + N - 1),u(t + N - 1)) =
P1(t,... ,t + N - 1)
P2(t,... ,t + N - 1)
= (Γ1(γ) Γ2(γ) . . . Γp(γ))
,
Pp(t,... ,t + N - 1)
где
Pi(t,... ,t + N - 1) := Pi(y(t),u(t),... ,y(t + N - 1),u(t + N - 1)).
26
Имея последнее представление в виду, запишем, как и в линейном случае [1],
выражение для значения величины y для T моментов времени:
(2.2)
YT (t) = (y(t) y(t-1) ... y(t-T +1)) = (Γ1(γ) Γ2(γ) ... Γp(γ))PT
(t),
P1(t - 1, . . . , t - N) P1(t - 2, . . . , t - N - 1) . . . P1(t - T, . . . , t - N - T + 1)
P2(t - 1, . . . , t - N) P2(t - 2, . . . , t - N - 1) . . . P2(t - T, . . . , t - N - T + 1)
PT (t) =
.
Pp(t - 1, . . . , t - N) Pp(t - 2, . . . , t - N - 1) . . . Pp(t - T, . . . , t - N - T + 1)
Характерной особенностью полученного выражения является то, что па-
раметры, значения которых могут быть неизвестны, находятся в строке
1(γ) Γ2(γ) . . . Γp(γ)), строка YT (t) и матрица PT (t) зависят только от
измеряемых значений векторов управления и выхода и элементов вектора γ,
значения которых известны. Последнее позволяет осуществлять диагности-
рование с использованием строки YT (t) и матрицы PT (t), т.е. без знания зна-
чений ряда или всех параметров системы.
2.2. Генерация невязки
Для генерации невязки, на основе которой принимается решение о возник-
ших дефектах, используется несколько методов [1], остановимся на методе,
основанном на анализе ядра матрицы PT (t). В [2, 5] размер временного окна T
выбирается минимальным, при котором rank (PT (t)) = rank (PT-1(t)). Отсю-
да следует, что последний столбец матрицы PT (t) линейно выражается через
предыдущие столбцы, т.е. существует ненулевой вектор v(T ), для которого
PT (t)v(T) = 0; это означает, что матрица PT (t) имеет непустое ядро. Из (2.2)
тогда следует YT (t)v(T ) = 0, которое является соотношением избыточности.
Из сказанного следует, что правило генерации невязки
(2.3)
rT (t) = YT (t)v(T), v(T) ker(PT
(t)),
рассмотренное для непрерывных систем в [2, 5], является робастным в том
смысле, что оно строится без знания значений некоторых параметров систе-
мы. Вычисление элементов соотношения (2.3), однако должно производиться
непосредственно в процессе диагностирования в реальном масштабе време-
ни. Этим правило (2.3) отличается от известных методов генерации невязки,
опирающихся на полное описание системы, что дает возможность выполнить
ряд вычислений предварительно - перед проведением процесса диагностиро-
вания. Отметим, что невязка rT (t) генерируется для каждого момента t, т.е.
временное окно является скользящим.
Недостатком рассмотренного подхода является то, что при наличии по-
грешностей в результатах измерений равенство невязки rT (t) нулю не будет
выполняться даже при отсутствии дефектов. Этот недостаток можно устра-
нить за счет использования пороговой логики принятия решений путем срав-
нения невязки с неким адаптивным порогом, некоторые правила построения
которого приведены в [8].
Ясно, что операции определения T , при котором rank (PT (t)) =
= rank(PT-1(t)), требуют определенных ресурсов времени и не всегда мо-
27
гут быть реализованы в реальном масштабе времени. Кроме того, вычисле-
ние рангов матриц, сформированных на основе экспериментальных данных,
является плохой обусловленной задачей. Для устранения этих недостатков
можно выбирать T заранее так, чтобы число столбцов матрицы PT (t) бы-
ло больше числа ее строк, так будет при T p + 1; в этом случае заведомо
выполняется равенство rank (PT (t)) = rank (PT-1(t)).
2.3. Поиск дефектов
При решении задачи поиска дефектов необходимо построить банк преоб-
разованных моделей, каждая из которых будет чувствительна к одной груп-
пе дефектов и нечувствительна к остальным. Детальный анализ дефектов
с этой точки зрения содержится в [8], не останавливаясь на нем, уточним,
что модель, нечувствительная к дефекту ρi, строится на основе результатов
рассматриваемой ниже теоремы 2.
Отметим, что модель, нечувствительная к дефекту ρi, может оказаться
нечувствительной и к некоторым другим дефектам, это можно выяснить по
виду системы уравнений (1.3) или (1.4): если в ее описание не входят па-
раметры γi, γj , . . . , γk, то модель нечувствительна к дефектам ρi, ρj , . . . , ρk и
чувствительна к остальным. Для отражения этих особенностей и принятия
решения о возникшем дефекте используется матрица синдромов S. Для па-
раметра γi строится вектор-синдром Si по следующему правилу: если этот
параметр входит в описание k-й модели, полагаем Ski = 1, в противном слу-
чае Ski = 0 (см. пример).
Совпадение синдромов Si и Sj для некоторых параметров γi и γj означает,
что дефекты ρi и ρj в рамках рассматриваемого подхода неразличимы. В этом
случае параметры γi и γj помещаются в одно множество Kd для некоторо-
го d. Параметры, имеющие индивидуальные синдромы, образуют множества,
включающие в себя только один параметр. В результате компоненты вектора
параметров (γ1, γ2, . . . , γs)T раскладываются в семейство непересекающихся
множеств K = {K1, K2, . . . , Kw}. Матрица синдромов S строится из множе-
ства векторов-синдромов, соответствующих элементам множества K, как из
столбцов.
3. Преобразование исходной модели
Условия и способы приведения исходной системы к виду (1.3) или (1.4),
рассматриваемые ниже, опираются на математический аппарат алгебры
функций, который детально изложен в [7]. Основные конструкции этого ап-
парата - отношение частичного порядка на множестве функций, бинарные
операции × и, бинарное отношение Δ, а также операторы M и m - кратко
описаны в Приложении.
3.1. Общий результат
Теорема 1. Система (1.1) преобразуема к виду (1.3) или (1.4) тогда и
только тогда, когда
(3.1)
m(h × (m(h × . . . × m(h) . . .))) ⊕ h = 1
(оператор m используется N раз).
28
Доказательство. Необходимость. Пусть ϕ(1),...,ϕ(N) и ψ - функ-
ции с необходимыми свойствами, т.е. преобразованные координаты (1.2) под-
чиняются уравнениям (1.3) (для конкретности рассмотрим только эти урав-
нения). Заменим в последнем из них переменную x∗N (t + 1) на ϕ(N)(x(t + 1))
и y(t) на h(x(t)):
ϕ(N)(x(t + 1)) = ϕ(N)(f(x(t), y(t))) =
(h(x(t)), u(t)).
По определению отношения Δ и оператора m отсюда следует включение
(h, ϕ(N)) Δ и неравенство ϕ(N) m(h). По аналогии из второго уравнения
в (1.3) следует неравенство
(
)
(3.2)
ϕ(i) m h × ϕ(i+1) × . . . × ϕ(N)
,
i = N - 1,...,2.
C учетом свойств операции × и оператора m получаем
при i = N - 1
(
)
ϕ(N-1) m h × ϕ(N)
m(h × m(h));
при i = N - 2
(
)
ϕ(N-2) m h × ϕ(N-1) × ϕ(N)
m(h × m(h) × m(h × m(h))).
Так как
h h × m(h),
то
m(h) m(h × m(h)),
а тогда
m(h) × m(h × m(h))= m(h × m(h)),
что в итоге дает
ϕ(N-2) m(h × m(h × m(h))).
Продолжая аналогично, получим в итоге из (3.2)
ϕ(1) m(h × (m(h × . . . × m(h) . . .)))
(оператор m используется N раз). Запишем равенство y = ψ(y) при
y = h(x1) = h(ϕ(1)(x)) и y = h(x) в виде h(ϕ(1)) = ψ(h). Поскольку h и
ψ - нетривиальные функции, то из последнего равенства имеем ϕ(1) ⊕ h = 1,
откуда с учетом полученного выше выражения для функции ϕ(1) следует
неравенство (3.1).
29
Достаточность. На основе неравенства (3.1) можно задать различные
функции ϕ(1), . . . , ϕ(N) и, следовательно, получить различные модели без об-
ратных связей. Для простоты рассмотрим случай приведения исходной мо-
дели к (1.4). Положим
ϕ(N) = m(h), x∗N = ϕ(N)(x),
ϕ(i) = m(h × ϕ(i+1)), x∗i = ϕ(i)(x), i = N - 1, . . . , 1.
Заметим, что если при некотором i функция m(h × ϕ(i+1)) является вектор-
ной, то для уменьшения размерности строящейся модели в качестве ϕ(i)
можно выбрать скалярную функцию такую, что ϕ(i) m(h × ϕ(i+1)); ска-
занное относится и к выбору функции ϕ(N). Отметим, что такой выбор не
всегда приводит к требуемому результату, в этом случае необходимо выби-
рать векторную функцию m(h × ϕ(i+1)). По определению оператора m имеем
(α, m(α)) Δ для произвольной функции α, поэтому
(
)
(h, m(h)) Δ и h × ϕ(i+1), m(h × ϕ(i+1))
Δ,
что после подстановки дает
(
)
(
)
h, m(ϕ(N))
Δ и h×ϕ(i+1)(i)
Δ, i = N - 1,... ,1.
По определению отношения Δ это означает, что для некоторых функций
,...,
справедливы уравнения (1.4). Из определения функций ϕ(1), . . .
...,ϕ(N) получаем
ϕ(1) = m(h × (m(h × . . . × m(h) . . .)))
(оператор m используется N раз). Отсюда и из неравенства (3.3) следует
существование функций h и ψ таких, что
h(ϕ(1)) = ψ(h) = ϕ(1) ⊕ h = 1,
или
y = h(ϕ(1)(x)) = ψ(y).
Отметим, что теорема справедлива и для систем с непрерывным временем,
это следует из результатов [7].
Пример 1. Рассмотрим дискретную систему
x1(t + 1) = (x2(t)(x3(t) + x5(t)) + x4(t))x4(t)u(t)),
x2(t + 1) = x1(t)/x4(t),
x3(t + 1) = x1(t)x5(t)/x4(t) - |x3(t) + x5(t)| - u(t),
x4(t + 1) = x2(t)(x3(t) + x5(t)) + x4(t),
x5(t + 1) = |x3(t) + x5(t)| + u(t),
y1(t) = x5(t), y2(t) = x4(t).
30
Проверим выполнение условия (3.1) для h(x) = (x5, x4)T : так как
(h × u) ⊕ f = x4 = f1/f4,
то
m(h) = x1/x4.
Аналогичные вычисления дают следующее:
m(h × m(h)) = (x1/x4, x2, x3 + x5)T,
m(h × m(h × m(h))) = (x1/x4, x2, x3 + x5, x4, x5)T.
Нетрудно видеть, что при N = 3 условие (3.1) выполняется, при этом
m(h × m(h × m(h))) h.
Положим
ϕ(3)(x) = m(h) = x1/x4.
Так как
m(h × x1/x4) = (x1/x4, x2, x3 + x5)T,
примем
ϕ(2)(x) = (x2, x3 + x5)T,
а поскольку
m(h × ϕ(2)) = (x1/x4, x2, x3 + x5, x4, x5)T,
положим
ϕ(1)(x) = (x4, x5)T.
Отметим, что выбор в качестве ϕ(i) скалярных функций приводит к невоз-
можности построения преобразованной модели, о чем было сказано в дока-
зательстве достаточности.
Положим x5 = x1/x4, x4 = x2, x3 = x3 + x5, x2 = x5, x1 = x4. В резуль-
тате преобразованная модель принимает вид
x1(t + 1) = x3(t)x4(t) + y2(t),
x2(t + 1) = |x3(t)| + u(t),
x3(t + 1) = x5(t)y1(t),
x4(t + 1) = x5(t),
x5(t + 1) = y2(t)u(t),
y1(t) = x2(t), y2(t) = x1(t).
31
Проведя серию временных сдвигов и подстановок на основе полученной
модели, нетрудно получить ее вход-выходное описание:
y1(t + 3) = |y1(t + 1)y2(t)u(t)| + u(t + 2),
y2(t + 3) = y2(t + 2) + y1(t + 1)(y2(t)u(t))2.
Отметим, что получение аналогичного вход-выходного описания путем вре-
менных сдвигов и подстановок на основе исходной модели - нетривиальная
задача из-за громоздкости получаемых выражений и необходимости их де-
тального анализа. Так, например, заменяя x5(t) на y1(t) и x4(t) на y2(t), по-
лучаем
y2(t + 1) = x2(t)(x3(t) + y1(t)) + y2(t),
)
x1(t)
(x1(t)y1(t)
y2(t + 2) =
- |x3(t) + y1(t)| - u(t) + y1(t + 1)
+ y2(t + 1);
y2(t)
y2(t)
ясно, что только сложный анализ последнего выражения даст нужный ре-
зультат.
3.2. Преобразование к модели, свободной от параметра
Для решения задачи локализации дефектов введем для i-й компоненты
вектора параметров γ векторную функцию α(i) с максимальным числом ком-
понент такую, что композиция α(i)(f(x, u, γ)) не содержит параметр γi. Если
функция f дифференцируема по этому параметру, сформулированное усло-
вие может быть записано в виде
α(i)(f(x,u,γ)) = 0, i = 1,... ,s.
∂γi
Известно [8], что модель, заданная в новых координатах с помощью преоб-
разования ϕ = ϕ(1) × ϕ(2) × . . . × ϕ(N), не будет зависеть от i-й компоненты
вектора параметров, т.е. будет нечувствительной к i-му дефекту, если вы-
полняется следующее условие:
(3.3)
α(i)
ϕ.
Условия построения модели, заданной в новых координатах и не зависящей
от i-й компоненты вектора параметров, даются следующей теоремой.
Теорема 2. Система (1.1) преобразуема к модели (1.4) c тождествен-
ной функцией h, не зависящей от i-й компоненты вектора параметров,
тогда и только тогда, когда выполняются следующие условия:
(3.4)
α(i)
m(h) = 1,
(
(
(
(
)
)))
(
)
(3.5)
m h × α(i)m h × ... × α(i)m(h) ...
⊕ α(i) ⊕h
=1
(оператор m используется N раз).
32
Доказательство. Необходимость. Пусть система (1.1) преобра-
зуема к модели (1.4), не зависящей от i-й компоненты вектора парамет-
ров. Это означает, что существуют функции ϕ(1), ϕ(2), . . . , ϕ(N), ψ, реализую-
щие преобразования (1.2), при этом α(i) ϕ = ϕ(1) × ϕ(2) × . . . × ϕ(N). Из по-
следнего условия получаем α(i) ϕ(j) для всех j. Из теоремы 1 следует
(h × ϕ(j+1), ϕ(j))) Δ, j = 1, . . . , N - 1, (h, ϕ(N)) Δ, что с учетом свойств
оператора m и неравенства α(i) ϕ(j) дает ϕ(j) α(i) m(h × ϕ(j+1)) и
ϕ(N) α(i) m(h) соответственно. Из последнего неравенства следует усло-
вие (3.4). Рассмотрим для простоты случай N = 3, тогда
(
)
(
(
))
ϕ(2) α(i) m h × ϕ(3)
α(i)m h × α(i)m(h)
,
(
)
(
(
(
(
))))
(3.6)
ϕ(1) α(i) m h × ϕ(2)
α(i) m h× α(i)m h× α(i) m(h)
Так как y = h(x1) и y = φ(y), то (φ(h))(x) = (h(ϕ(1)))(x). Это означает,
что h ⊕ ϕ(1) = 1, а тогда с учетом (3.6) получаем условие (3.5).
Достаточность. Из (3.5) следует существование функций γ1, γ2 и ϕ(1)
таких, что
(
(
(
(
(
)
))))
ϕ(1) = γ1 m h × α(i) m h × . . . × α(i) m(h) . . .
=
(
)
=γ2
α(i) ⊕ h
= 1.
Из равенства ϕ(1) = γ2(α(i) ⊕ h) следует ϕ(1) h и ϕ(1) α(i). Из определения
функции ϕ(1), условия (3.5) и свойств операторов m и M получаем
(
(
(
)
))
(
)
h × α(i)m h × ... × α(i)m(h) ...
M ϕ(1)
Выберем в качестве ϕ(2) максимальную функцию, удовлетворяющую усло-
виям h × ϕ(2) M(ϕ(1)) и
(
(
)
)
α(i) m h × ... × α(i) m(h) ...
ϕ(2)
(оператор m используется N - 1 раз), тогда справедливо включение
(h × ϕ(2), ϕ(1)) Δ. Из предыдущего следует, что это можно сделать всегда,
положив
(
(
)
)
ϕ(2) = α(i) m h × . . . × α(i) m(h)
Тогда
(
(
)
)
ϕ(2) m h × . . . × α(i) m(h) . . . и ϕ(2) α(i).
По аналогии строятся функции ϕ(3), . . . , ϕ(N-1), для которых справедливы
включение (h × ϕ(j+1), ϕ(j)) Δ и неравенство ϕ(j) α(i), j = 3, . . . , N - 1.
33
Используя оператор M N - 1 раз, получим из (3.5) нетривиальную (в си-
лу условия (3.4)) функцию ϕ(N) = α(i) m(h), откуда следует ϕ(N) m(h)
и ϕ(N)α(i), т.е. справедливо включение (h,ϕ(N)) Δ. Так как ϕ = ϕ(1)×
×ϕ(2) × ... × ϕ(N) и ϕ(j) α(i) для всех j, то ϕ α(i), т.е. справедливо усло-
вие (3.3). Из всего сказанного и определения отношения Δ следует, что для
некоторых функций
,...,
справедливы уравнения (1.4). Из неравен-
ства ϕ(1) h следует существование функции ψ такой, что ϕ(1) = ψ(h), отку-
да y = ψ(y) = ψ(h(x)) = ϕ(1)(x), т.е. h - тождественная функция.
Пример 2. Рассмотрим дискретную систему
x1(t + 1) = γ2x2(t) + x3(t) - x2(t),
x2(t + 1) = u(t)|x1(t)|,
x3(t + 1) = u(t)|x1(t)| + x1(t)x4(t),
x4(t + 1) = γ1x5(t)(x3(t) - x2(t)),
x5(t + 1) = u(t)x1(t),
y1(t) = x1(t), y2(t) = x3(t) - x2(t).
Проверим выполнение условий теоремы 2 для первого дефекта, которо-
му соответствует функция α(1) = (x1, x2, x3, x5)T. Так как m(h) = (x2, x5)T,
то α(1) m(h) = (x1, x2, x3, x5)T (x2, x5)T = (x2, x5)T = 1, т.е. выполняется
условие (3.4). Далее
(
)
h × α(1)m(h)
= (x1, x2, x3, x5)T
и
(
)
m h × (α(1)m(h))
= (x1, x2, x4, x5)T,
а так как α(1) ⊕ h = (x1, x3 - x2)T, то выполняется условие (3.5).
Положим ϕ(1)(x) = (x1, x2, x4, x5)T (x1, x3 - x2)T = x1; так как ϕ(2)
должно удовлетворять условию
ϕ(2) α(1) m(h) = (x1, x2, x3, x5)T (x2, x5)T = (x2, x5)T,
примем ϕ(2)(x) = x2.
Положим x1 = x1, x2 = x2; преобразованная модель, нечувствительная к
первому дефекту, имеет вид
x1(t + 1) = γ2x2(t) + y2(t),
x2(t + 1) = u(t)|y1(t)|,
y1(t) = x1(t).
Для второго дефекта, которому соответствует функция
α(2)(x) = (x2,x3,x4,x5)T,
34
аналогичным образом можно показать, что x1 = x2, x2 = x3, x3 = x4,
x4 = x5; модель, нечувствительная к этому дефекту, принимает вид
x1(t + 1) = u(t)|y1(t)|,
x2(t + 1) = u(t)|y1(t)| + y1(t)x3(t),
x3(t + 1) = γ1x4(t)y2(t),
x4(t + 1) = u(t)y1(t),
y2(t) = x2(t) - x1(t).
На основе этих моделей нетрудно получить вход-выходное описание системы.
4. Решение на основе логико-динамического подхода
Теоремы 1 и 2 дают исчерпывающий ответ на вопрос о возможности пред-
ставления исходной системы в виде, необходимом для применения непара-
метрического метода. Они дают общее решение задачи, когда преобразова-
ние (1.2) исходной системы к форме (1.3) или (1.4) может быть нелинейным,
и требуют применения сложного математического аппарата алгебры функ-
ций. В случае ограничения класса преобразований линейными функциями
решение может быть получено методами линейной алгебры на основе логико-
динамического подхода [8]. Для возможности его применения исходная систе-
ма должна быть представлена в виде
ψ1(A1x(t),u(t))
x(t + 1) = F x(t) + Gu(t) + C
+Didi(t),
(4.1)
ψq(Aqx(t),u(t))
i=1
y(t) = Hx(t),
где x ∈ Rn, u ∈ Rm, y ∈ Rl - векторы состояния, управления и выхода;
F, G и C - матрицы при номинальных значениях параметров, F и G описы-
вают линейную динамику; H, D1, . . . , Ds - известные постоянные матрицы,
d1(t),... ,ds(t) - скалярные функции, описывающие дефекты: если i-й па-
раметр имеет номинальное значение, то di(t) = 0, в противном случае di(t)
становится неизвестной функцией времени, i = 1, . . . , s; ψ1, . . . , ψq - произ-
вольные нелинейности, A1, . . . , Aq - матрицы-строки. Модель (4.1) всегда мо-
жет быть получена из нелинейной модели (1.1) путем ряда эквивалентных
преобразований, детально описанных в [8].
Известно [8], что логико-динамический подход реализуется в три этапа,
на первом из которых из модели (4.1) удаляется нелинейная составляющая,
на втором рассматриваемая задача решается для линейной системы с до-
полнительным ограничением линейного же характера, на третьем к полу-
ченному линейному решению добавляется преобразованная нелинейная со-
ставляющая. Коротко опишем решение на втором этапе, в основе которого
лежит процедура синтеза модели, представленной в канонической форме [8],
на основе которой нетрудно получить вход-выходные соотношения, лежащие
в основе непараметрического метода.
35
Приведем описание преобразованной модели:
x∗i(t + 1) = x∗i+1(t) + Jiy(t) + G∗iu(t), i = 1,... ,k - 1,
(4.2)
x∗k(t + 1) = Jky(t) + G∗ku(t),
y(t) = x1(t),
где x Rk - вектор состояния модели, x = Φx, y = Ry. Известно, что мат-
рицы, описывающие исходную систему и эту модель, в также матрицы Φ и R
удовлетворяют следующим уравнениям [8]:
(4.3)
RH = Φ1, ΦiF = Φi+1 + JiH, i = 1,... ,k - 1, ΦkF = Jk
H,
где Φi и Ji - i-е строки матриц Φ и J, i = 1, . . . , k, k - размерность модели.
Дополнительное ограничение на матрицу Φ имеет вид
(
)
Φ
(4.4)
A=A
H
для некоторой матрицы A, где
(
A=
AT1 ... ATq
)T.
Равенство (4.4) эквивалентно условию [8]
(
)
(
)
(4.5)
rank
ΦT HT
= rank
ΦT HT AT
Известно, что уравнения (4.3) можно свернуть в одно:
(4.6)
RHFk = J1HFk-1 + J2HFk-2 + ... + Jk
H.
Решение этого уравнения дает минимальную размерность k и матрицы R и J.
Для построения модели, не содержащей i-й компоненты вектора парамет-
ров, напомним, что соответствующее условие имеет вид ΦDi = 0 [8]. Учет
этого условия приводит к уравнению [1]
(4.7)
(R
-J1
-J2
-Jk) (U(k) B(k)
) = 0,
где
HFk
HFk-1
U(k) =
,
HFk-2
H
(4.8)
HDi HFDi HF2Di ... HFk-1D
i
0
HDi HFDi ... HFk-2Di
B(k) =
.
0
0
HDi ... HFk-3Di
0
0
0
0
36
Минимальное значение размерности k, при котором это уравнение имеет
нетривиальные решения, определяется условием
(4.9)
rank (U(k) B(k)
) < l(k + 1).
При его выполнении найдется вектор-строка (R
-J1
-J2
-Jk) , удо-
влетворяющая уравнению (4.7). Далее по формулам (4.3) определяются стро-
ки матрицы Φ и проверяется условие (4.5), при его выполнении рассчитыва-
ется матрица G = ΦG и строится модель в форме (4.2). Если условие (4.5) не
выполняется, рекомендуется найти другое решение уравнения (4.7) или уве-
личить размерность k; невыполнение этого условия при всех k < n означает,
что модели, не содержащей i-й компоненты вектора параметров, не существу-
ет. Далее будем полагать, что условие (4.5) выполняется, т.е. уравнение (4.4)
имеет решение.
Для построения нелинейной составляющей заметим, что она описывается
выражением
ψ1(A1z(t),u(t))
(4.10)
C
,
ψq(A∗qz(t),u(t))
(
)
(
)
x
Φ
где z =
, строка A1 определяется из уравнения Ai = A∗i
,
y
H
i = 1,...,q, соответствующего (4.4), C = ΦC. Если выражение (4.10) имеет
вид, соответствующий правым частям моделей (1.3) или (1.4), т.е. не содер-
жит обратных связей, составляющая (4.10) добавляется к ранее построенной
линейной части. В противном случае следует найти другое решение уравне-
ния (4.7) при прежней или увеличенной размерности.
Отметим, что если все строки матрицы A линейно выражаются через стро-
ки матрицы H, то выражение (4.10) заведомо имеет требуемый вид. Если это
условие не выполняется, добиться его выполнения можно путем расширения
вектора выхода исходной системы. Детальное рассмотрение этого вопроса,
однако, выходит за рамки настоящей работы.
Дальнейшие действия с полученной моделью не отличаются от тех, кото-
рые были описаны в разделе 2.
В некоторых случаях система не может быть преобразована к виду (1.3)
или (1.4), однако с помощью эвристических приемов она может быть приве-
дена к виду (2.1), рассмотрим это на примере.
Пример 3. Пусть преобразованная модель описывается системой урав-
нений
x1(t + 1) = x2(t) + η1(x1(t)) + J1y(t),
x2(t + 1) = x3(t) + η2(x2(t)) + J2y(t),
x3(t + 1) = η3(x3(t)) + J3y(t),
y(t) = x1(t),
37
где η1 ÷ η3 — произвольные нелинейные функции, вектор управления для
простоты опущен, элементы матриц-строк J1, J2 и J3 зависят от парамет-
ров исходной системы. Особенность модели состоит в том, что правая часть
уравнения для каждой переменной содержит нелинейную функцию той же
переменной. Рассмотрим ряд временных сдвигов и подстановок для перемен-
ной y:
y(t+2) = x1(t + 2) = x2(t + 1) + η1(x1(t + 1)) + J1y(t + 1) =
= η1(y(t + 1)) + J1y(t + 1) + x3(t) + η2(x2(t)) + J2y(t),
y(t+3) = η1(y(t+2))+J1y(t+2)+x3(t+1)+η2(x2(t+1))+J2y(t+1) =
= η1(y(t + 2)) + J1y(t + 2) + J2y(t + 1)+
+ η2(x3(t) + η2(x2(t)) + J2y(t)) + η3(x3(t)) + J3y(t).
Из выражения для y(t + 2) найдем сумму x3(t) + η2(x2(t)) + J2y(t):
x3(t) + η2(x2(t)) + J2y(t) = y(t + 2) - η1(y(t + 1)) - J1y(t + 1)
и подставим ее в уравнение для y(t + 3) вместо аргумента функции η2:
y(t + 3) = η1(y(t + 2)) + J1y(t + 2) + J2y(t + 1)+
+ η2(y(t + 2) - η1(y(t + 1)) - J1y(t + 1)) + η3(x3(t)) + J3y(t).
Если далее из последнего соотношения выразить сумму η3(x3(t)) + J3y(t),
найти выражение для y(t + 4) и подставить эту сумму вместо аргумента
функции η3, то итоговое выражение для y(t + 4) будет содержать только
временные сдвиги векторов y и y, что и требуется для применения непара-
метрического подхода.
5. Практический пример
Рассмотрим модель электропривода при отсутствии внешнего нагрузочно-
го момента и без учета вязкого трения:
x1 =
1x2,
ip
KM
x2 =
x3 -
Kd sign(x2),
JH
JH
Kω
R
KU
x3 = -
x2 -
x3 +
u,
L
L
L
где x1 - угол поворота выходного вала редуктора, x2 - угловая скорость вра-
щения вала электродвигателя, x3 - ток электродвигателя, ip - передаточное
отношение редуктора, JH - момент инерции ротора двигателя и вращающих-
ся частей редуктора, приведенный к ротору, KM - моментный коэффициент
38
электродвигателя, Kω - коэффициент противо-э.д.с., Kd - момент сухого тре-
ния в двигателе, R - активное сопротивление обмотки якоря электродвига-
теля, L - индуктивность якорной обмотки электродвигателя, KU - коэффи-
циент усиления усилителя мощности.
Дискретизированную модель представим в виде
x+1 = γ1x2 + x1,
(5.1)
x+2 = γ2x3 + γ3sign(x2) + x2,
x+3 = γ4x2 + γ5x3 + γ6u,
где для простоты используются обозначения x+ = x(t + 1), x = x(t), u = u(t),
коэффициенты γ1 ÷ γ6 представляют описанные выше параметры и интервал
дискретизации.
Из уравнений (5.1) ясно, что коэффициенты γ2, γ3, а также γ4 ÷ γ6 одина-
ковым образом влияют на переменные x2 и x3 соответственно, поэтому можно
говорить о трех дефектах. Полагая, что измеряемыми являются переменные
x1 и x3, рассматриваемую модель опишем следующими матрицами:
(
)
1
γ1
0
0
1
0
0
F = 0
1
γ2
,G=0
,H=
,
0
0
1
0
γ4
γ5
γ6
0
ϕ(x, u) = sign(Ax), A = (0 1 0), C = γ3
,
0
1
0
0
D1 = 0
,D2 = 1
,D3 = 0
.
0
0
1
Следуя логико-динамическому подходу, удалим нелинейный член из мо-
дели и рассмотрим ее линейную часть. Непосредственные вычисления для
первого дефекта дают следующее:
(
)
1
HD1 =HFD1 =
,
0
1
2γ1
γ1γ2
1
1
0
γ4(1 + γ5) γ2γ4 + γ24
0
0
1
γ1
0
0
1
(U(2), B(2)) =
0
γ4
γ5
0
0
1
0
0
0
0
0
0
1
0
0
Ясно, что rank (U(2), B(2)) < 2(2 + 1) = 6, поэтому уравнение (4.7) имеет ре-
шение, приведем его в виде
(R
-J1
-J2) = (0
1
0
-(1 + γ5)
0
-(γ2γ4 - γ5)) .
39
Тогда
(
)
γ6
Φ1 = (0 0 1), Φ2 = (0 γ4
- 1), G =
6
Нетрудно проверить, что условие (4.5) выполняется и матрицу A можно
принять в виде A = (0 14 0 1), кроме того,
(
)
0
C = ΦC =
γ3γ4
В результате нелинейная преобразованная модель описывается уравне-
ниями
x+1 = x2 + (1 + γ5)y2 + γ6u,
x+2 = (γ2γ4 - γ5)y2 - γ6u + γ3γ4sign((y2 + x2)4),
y1 = x1 = y2,
где
x1 = x3, x2 = γ2x2 - x3.
Временные сдвиги дают следующее:
y2(t + 1) = x2(t) + (1 + γ5)y2(t) + γ6u(t),
y2(t + 2) = (γ2γ4 - γ5)y2(t) - γ6u(t) + γ3γ4sign((y2(t) + x2(t))4) +
+ (1 + γ5)y2(t + 1) + γ6u(t + 1).
Из выражения для y2(t + 2) следует, что сумма первых трех его слагаемых
из правой части может быть представлена как
y2(t + 2) - (1 + γ5)y2(t + 1) - γ6u(t + 1).
Имея это в виду, приведем выражение для y2(t + 3):
y2(t + 3) = (γ2γ4 - γ5)y2(t + 1) - γ6u(t + 1) + (1 + γ5)y2(t + 2) + γ6u(t + 2) +
+ γ3γ4sign((y2(t + 1) + y2(t + 2) - (1 + γ5)y2(t + 1) - γ6u(t + 1))4).
Поскольку функция sign содержит параметры γ4, γ5 и γ6, то для применения
непараметрического метода их значения должны быть известны.
Выражения для остальных сдвигов могут быть получены по аналогии. Так
как последнее выражение содержит 5 слагаемых, то T = 6. Обозначим функ-
цию sign() в последнем выражении через z(t) и запишем его в виде (3.2),
приняв y = Ry = y2 :
Y6(t) = (y2(t) y2(t - 1) y2(t - 2) y2(t - 3) y2(t - 4) y2(t - 5)) =
= (γ2γ4 - γ5
6
1+γ5
γ6
γ3γ4) P6(t),
40
r
0,14
0,12
0,10
0,08
0,06
0,04
0,02
0
t
0
20
40
60
80
100
Рис. 1. Поведение невязки r(1)6
при изменении параметров k1 и k2.
r
0,35
0,30
0,25
0.20
0,15
0,10
0,05
t
0
0
20
40
60
80
100
Рис. 2. Поведение невязки r(1)6
при изменении параметров k3 и k4.
41
где
y2(t - 2) y2(t - 3) y2(t - 4) y2(t - 5) y2(t - 6) y2(t - 7)
u(t - 2) u(t - 3) u(t - 4) u(t - 5) u(t - 6) u(t - 7)
P6(t) =
y2(t - 1) y2(t - 2) y2(t - 3) y2(t - 4) y2(t - 5) y2(t - 6)
.
u(t - 1) u(t - 2) u(t - 3) u(t - 4) u(t - 5) u(t - 6)
z(t - 3) z(t - 4) z(t - 5) z(t - 6) z(t - 7) z(t - 8)
Невязка формируется в виде r(1)6(t) = Y(1)6(t)v(6), где v(6) ker(P(1)6(t)).
Можно показать, что уравнение (4.7) для второго дефекта не имеет реше-
ний. Выполняя требуемые действия для третьего дефекта, получаем следую-
щую преобразованную модель:
x+1 = x2 + 2y1,
x+2 = -y1 + γ1γ2y2 + γ1γ3sign((y1 + x2)1),
y1 = x1 = y1,
где x1 = x1, x2 = γ1x2 - x1. Произведя временные сдвиги, можно заметить,
что в результирующем выражении функция sign содержит параметр γ1, по-
этому для применения непараметрического метода его значение должно быть
известно. Невязка r(2) формируется по аналогии с r(1)6.
Нетрудно видеть, что матрица синдромов в рассматриваемом примере
имеет вид
( 0 1 1)
S=
,
1
1
0
где строки соответствуют невязкам r(1)6 и r(2), столбцы - множествам
K1 =1}, K2 =22}, K3 =456}.
Для моделирования примем k1 = k2 = k6 = 1, k3 = k4 = k5 = -1, u(t) =
= 2sin(t/10). На рис. 1 представлены результаты моделирования в случае
k1 = 1,2 в момент t = 40 и k2 = 1,2 в момент t = 70; видно, что невязка r(1)6
нечувствительна к изменению первого параметра и чувствительна ко второ-
му. На рис. 2 представлены результаты моделирования в случае k3 = -1, 2 в
момент t = 40 и k4 = -1, 2 в момент t = 70; видно, что невязка чувствительна
к изменению обоих параметров.
6. Заключение
В работе рассмотрен подход к решению задачи диагностирования техни-
ческих систем, описываемых нелинейными моделями, путем проверки соот-
ношений, существующих между входами и выходами системы, измеряемыми
на конечном интервале времени. Получены критерии возможности построе-
ния таких соотношений для заданной нелинейной системы. Для проверки
используется правило, реализуемое на основе непараметрического метода.
Теоретические результаты иллюстированы примером диагностирования ре-
альной технической системы.
42
ПРИЛОЖЕНИЕ
Используемый в работе математический аппарат (алгебра функций) со-
держит четыре основные конструкции: 1) отношение частичного предпоряд-
ка, обозначаемое как, 2) две бинарные операции × и, 3) бинарное от-
ношение Δ, 4) операторы m и M. Первые два элемента определены на мно-
жестве VW произвольных векторных функций с некоторой областью опре-
деления W , в то время как два последних - на множестве VX векторных
функций с областью определения в пространстве состояний X. Рассмотрим
эти элементы более подробно.
1. Отношение частичного предпорядка. Для произвольных функций
α, β ∈ VW записывается α β, если существует функция γ такая, что
β(w) = γ(α(w))
∀w ∈ W. Определение означает, что каждая компонента
функции β может быть выражена через компоненты функции α. Если α β
и βα, то функции α и β называются эквивалентными, что обозначается
как α= β.
Отметим, что= является отношением эквивалентности. Оно разбива-
ет множество VW на классы эквивалентности, содержащие эквивалентные
функции. Обозначим через VE множество классов эквивалентности, тогда от-
ношение является частичным порядком на этом множестве и пара (VE ,)
представляет собой решетку [9]. Существуют две специальные функции 0
и 1 такие, что для произвольной функции α выполняется условие 0 α 1.
Функция 0 эквивалентна тождественной функции, 1 - постоянной функции.
2. Бинарные операции. Поскольку пара (VE ,) представляет собой решет-
ку, то для произвольных функций α, β ∈ VW существуют две бинарные опе-
рации × and, определяемые следующим образом [9]:
α × β = inf(α,β), α ⊕ β = sup(α,β).
Указанные операции определяют функции с точностью до эквивалентности.
Правило вычисления операции × выглядит просто:
(
)
α(w)
(α × β)(w) =
β(w)
Операция может быть вычислена с использованием средств дифференци-
альной геометрии [7]. В простейших случаях может быть использовано опре-
деление операции α ⊕ β как наименьшей верхней грани функций α и β.
Рассмотрим пример вычисления функций α × β и α ⊕ β. Пусть W = R3,
)
(w1 +w2
(w1w3)
α(w) =
,
β(w) =
w3
w2w3
Тогда (α × β)(w)= (w1 + w2, w3, w1w3)T, (α ⊕ β)(w) = w3(w1 + w2).
3. Бинарное отношение Δ. Для заданных функций α, β ∈ VX
(α, β) Δ ⇐⇒ α(f(x, u, γ)) = f(β(x), u, γ)
43
для всех (x, u) ∈ X × U и некоторой функции f. Бинарное отношение Δ са-
мостоятельного значения не имеет и используется для определения операто-
ров.
4. Операторы m и M. Для произвольных функций α, β ∈ VX определяются
операторы m и M следующим образом.
Оператор m задает функцию m(α) ∈ SX , которая удовлетворяет усло-
виям: (i) (α, m(α)) Δ, (ii) если (α, β) Δ, то m(α) β.
Оператор M задает функцию M(β) ∈ SX , которая удовлетворяет усло-
виям: (i) (M(β), β) Δ, (ii) если (α, β) Δ, то α ≤ M(β).
Из последних определений следует, что для произвольной α выраже-
ние m(α) представляет собой минимальную функцию, формирующую пару
с α, и для произвольной β выражение M(β) представляет собой максималь-
ную функцию, формирующую пару с β.
Оператор m может быть вычислен следующим образом. Пусть функция γ
удовлетворяет условию (α × u) ⊕ f= γ(f). Тогда m(α)= γ. Можно пока-
зать, что функция γ всегда существует. Правило вычисления оператора M
не приводится, поскольку в работе оно не используется.
Операторы m и M обладают рядом свойств, которые детально описаны
в [7].
Отметим, что стандартные математические пакеты программ нельзя ис-
пользовать для вычисления операций и операторов алгебры функций. По-
этому на основе пакета Mathematica была разработана вычислительная сре-
да, доступная на сайте [10] и позволяющая организовать облачные вычисле-
ния для нахождения упомянутых операций и операторов, а также проверки
условий теорем 1 и 2. Характерной особенностью таких вычислений является
возможность их выполнения без наличия пакета Mathematica у пользователя
сайта.
СПИСОК ЛИТЕРАТУРЫ
1. Жирабок А.Н., Шумский А.Е., Павлов С.В. Диагностирование линейных дина-
мических систем непараметрическим методом // АиТ. 2017. № 7. С. 3-21.
Zhirabok A.N., Shumsky A.E., Pavlov S.V. Diagnosis of Linear Dynamic Systems
by the Nonparametric Method // Autom. Remote Control. 2017. V. 78. No. 7.
P. 1173-1188.
2. Ding S. Data-driven Design of Fault Diagnosis and Fault-tolerant Control Systems.
London: Springer-Verlag, 2014.
3. Lindner B., Auret L. Data-driven fault getection with process topology for fault
identification // Proc. 19th World IFAC Congr. 2014. P. 8903-8908.
4. Haghani A., Krueger M., Jeinsch T., Ding S., Engel P. Data-driven multimode
fault detection for wind energy conversion systems // Proc. 9th IFAC Sympos.
Safeprocesses’2015. P. 633-6385.
5. Шумский А. Функциональное диагностирование нелинейных динамических си-
стем с запаздыванием // АиТ. 2009. № 3. С. 172-184.
Shumskii A.E. Functional Diagnosis of Nonlinear Dynamic Delay Systems // Autom.
Remote Control. 2009. V. 70. No. 3. P. 502-512.
6. Shumsky A., Zhirabok A. Redundancy relations for fault diagnosis in hybrid
systems // Proc. 8th IFAC Sympos. Safeprocesses’2012. P. 1226-1231.
44
7. Жирабок А.Н., Шумский А.Е. Алгебраические методы анализа нелинейных ди-
намических систем. Владивосток: Дальнаука, 2008.
8. Шумский А.Е., Жирабок А.Н. Методы и алгоритмы диагностирования и от-
казоустойчивого управления динамическими системами. Владивосток: Изд-во
ДВГТУ, 2009.
9. Гретцер Г. Общая теория решеток. М.: Мир, 1982.
10. http://webmathematica.cc.ioc.ee/webmathematica/NLControl/main/index.html.
Статья представлена к публикации членом редколлегии С.А. Красновой.
Поступила в редакцию 12.02.2017
После доработки 15.04.2018
Принята к публикации 15.04.2018
45