CC BY
15
УДК 519.6
DOI: 10.14529/mmph190404
ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ КВАЗИЛИНЕЙНОГО РАЗНОСТНОГО УРАВНЕНИЯ
А.В. Панюков, Я.А. Мезал
Южно-Уральский государственный университет, г. Челябинск, Российская Федерация E-mail: paniukovav@susu.ru
Идентификация квазилинейного разностного уравнения сводится к задаче регрессионного анализа с взаимно зависимыми наблюдаемыми переменными. Это делает неэффективными классические схемы решения, основанные на методе наименьших квадратов и его вариациях. Нахождение оценок коэффициентов уравнения авторегрессии существенно осложняется плохой обусловленностью системы уравнений, представляющих собой необходимые условия минимума суммы квадратов отклонений. При этом оценки параметров задачи оказываются несостоятельными. Для решения подобных задач возможно применение обобщённого метода наименьших модулей (ОМНМ), сводимого к решению последовательности задач оценки коэффициентов уравнения регрессии по взвешенному методу наименьших модулей (ВМНМ). В статье предложен алгоритм решения задачи ВМНМ-оценивания, на основе ее сведения к задаче линейного программирования (ЛП) простой структуры. Простота структуры допустимого множества используемой задачи ЛП: пересечение линейного подпространства с параллелепипедом, - позволяет предложить эффективный алгоритм ее решения, основанный на методе проекции градиента. Алгебраическая вычислительная сложность предложенного алгоритма не превосходит величины O(N2M2), где N - количество коэффициентов в исследуемом уравнении, M -количество наблюдаемых значений. Данная оценка вычислительной сложности ВМНМ является наилучшей из известных.
Ключевые слова: метод наименьших модулей; модель авторегрессии; линейное программирование; метод проекции градиента; вычислительная сложность.
Введение
Проблема идентификации является в настоящее время одной из основных проблем системного анализа и его приложений в различных областях [1]. При идентификации предполагается экспериментальное изучение и сопоставление входных и выходных процессов, а задача идентификации состоит в выборе соответствующей математической модели. Модель должна быть такой, что ее реакция и реакция объекта на один и тот же входной сигнал должны быть, в известном смысле, близкими. Поиск модели, как правило, ведется в параметризованных классах, полученных «на основе физических закономерностей и представлений о процессах». Одним из таких классов является класс квазилинейных разностных уравнений. В качестве примеров можно привести:
• модели делового цикла [2]:
7, -Yt_i = ut + s{Yt_x - Jt_2) + v(Yt_l -It_2)3, где Yt - доход в момент времени t, ut - налоги (экзогенная переменная), v, s - параметры модели;
• модели оценки влияния различных факторов рыночной среды на процессы переноса (трансфера) и освоения инновационных (в том числе информационных) технологий [3]:
dy = u(t) + ay(t) + by2, dt
где y(t) - объём распространения инновации к моменту t; экзогенная функция u(t) и параметры а, b модели отражают суммарное число потенциальных потребителей инновационного продукта на рынке, степень внешних (влияние СМИ, рекламы, публикаций) и внутренних воздействий (непосредственное общение людей) на скорость адаптации;
• и др. [4]._
Панюков А.В., Мезал Я.А.
Параметрическая идентификация квазилинейного разностного уравнения
Как правило, все подобные модели имеют параметры, являющиеся коэффициентами систем интегро-дифференциальных уравнений. В данном случае исходная проблема идентификации может быть сведена к задаче идентификации коэффициентов разностного уравнения, т. е. к оценке их значений по наблюдаемым возмущенным траекториям. Такой подход правомерен при наличии большого объема априорной информации об искомых параметрах и при несоответствии реальных параметров принятой модели приводит к потере сходимости алгоритма.
В данной работе показано, что идентификация квазилинейного разностного уравнения сводится к задаче регрессионного анализа с взаимно зависимыми наблюдаемыми переменными. В первом разделе задача идентификации квазилинейного уравнения представлена в виде задачи ОМНМ-оценивания [5, 6], которая в соответствии с [7] может быть решена с помощью решения последовательности задач ВМНМ-оценивания, представляющих задачи кусочно-линейного программирования. В третьем разделе подробно рассмотрена задача ВМНМ-оценивания. Показано, что ее решение сводится к задаче линейного программирования (ЛП) простой структуры. Простота структуры допустимого множества используемой задачи ЛП: пересечение линейного подпространства с параллелепипедом, - позволяет предложить эффективный алгоритм ее решения, основанный на методе проекции градиента [8, 9]. Описание алгоритма и доказательство его результативности приведено в разделе 4.
Постановка задачи
Пусть е- значения переменных состояния (т.е. эндогенных переменных) в моменты времени I (еМ . Пусть {ха,хп,.. е М}^ - значения переменных управления (т.е.
экзогенных переменных) в моменты времени I (еМ . Пусть g.: <т —»К, 7=1,2,...т - заданные
функции, (е(
- случайные ошибки. В статье рассматривается проблема определения
неизвестных коэффициентов а^, а2, а3..., ат е М, I авторегрессии
т п
квазилинеинои модели
Г = \,2 ,...,М. (1)
м м
Введем новые обозначения с целью уменьшения громоздкости математических выражений
А =
а1 gl(Уt-к }т=1
а2 g2 {Уг-к }т=1
ат ь ; X = gm { Уг - к }т=1 х1г
ь2 х2 г
_ Ьп _ Хпг
г = 1,2,..., М; N = п + т.
В новых терминах задача (1) примет вид
у(=АТХ(+еп ¿ = 1,2,...,М. (2)
Наиболее известным методом определения коэффициентов уравнения регрессии является метод наименьших квадратов (МНК)
1 2
(3)
который является параметрическим методом, требующим выполнения ряда строгих ограничений: детерминированности переменных, независимости и нормальности распределения погрешностей измерений [10-12]. Даже незначительные нарушения этих предпосылок резко снижают эффективность оценок МНК [13, 14].
Если допустить погрешности в измеренных значениях эндогенных переменных у, г = 1,2,.,Т), то их наличие в значениях функций g .очевидно. Кроме того, эти
г=1
ошибки будут взаимно коррелированы и иметь распределения вероятностей, отличные от нормального распределения. Это делает неэффективными классические схемы решения, основанные на методе наименьших квадратов и его вариациях.
Нахождение оценок коэффициентов уравнения авторегрессии существенно осложняется плохой обусловленностью системы уравнений, представляющих собой необходимые условия минимума суммы квадратов отклонений, при этом оценки становятся несостоятельными. Альтернативой МНК является метод наименьших модулей (МНМ) [15, 16].
M
A = arg min УAT- yt I. (4)
ЛеЛ i=1
Возможными обобщениями МНМ являются взвешенный МНМ (ВМНМ)
м
^=argmin X(pt\ATXt-yt\\ pteR, t = \2,..M (5)
Aeh t=i
и обобщенный МНМ (ОМНМ)
м
Л' = arg min £p(|ЛТ- y,|), (6)
Лек i=1
где p(*) - выпуклая вверх дифференцируемая функция.
Основной проблемой при применении ВМНМ является отсутствие общих формальных правил выбора весовых коэффициентов. Идентификация квазилинейного уравнения (1) авторегрессии возможна применением ОМНМ с использованием установленной связи между ОМНМ- и ВМНМ-моделями [7], позволяющей решить проблему определения ОМНМ-оценок посредством итеративной процедуры ОМНМ-оценивания с ВМНМ-оценками [5, 6].
Проблема ВМНМ оценивания
Алгоритм ВМНМ-оценивания для задачи (5) приводит к решению задачи кусочно-линейного программирования
м
У р\ЛтX -yt\ ^ min, (7)
Luft | » st\ v /
i=i
для заданных Xt el", yt,Pt e R, t = \,2,...M . Задача (7) эквивалентна задаче линейного программирования
м
Zp.z. min , (8)
ATXt+zt>yt, t = l,2,...,M, (9)
-ЛтX + Ъ ^-У, t = !,2,—,M . (10)
Двойственной задаче (8)-(10) является задача
м
t=1
M
[и -V.
MXtj [и-vt) = 0, j = 1,2,...,N, (12)
и + V = Р, * = 1,2,..., М, (13)
иV >0, I = 1,2,...,М. (14)
Введем переменные = и -V, t = 1,2,...,М. Из условий (13)—(14) следует
р, + wt р, - wt
, , , = 1,2,.,М.
-р <^ < Р,, = 1,2,.,М. Поэтому оптимальное значение задачи (11)—(14) равно оптимальному значению задачи
М
тах, (15)
Панюков А.В., Параметрическая идентификация
Мезал Я.А. квазилинейного разностного уравнения
M
X Xwt = 0, j = 1,2,., N, (16)
t=i
-pt < w, < pt, t = 1,2,.,M. (17)
Если w* является оптимальным решением задачи (15)—(17), то оптимальное решение задачи (11)—(14) равно
. ft +w* . pt - w* ut =—-—, vt =—-—, t = 1,2,.,M.
Из условия комплементарности для пары взаимно двойственных задач (8)—(10) и (11)—(14) следует
(A*)T Xt = y„ t:-p, <w* <p„ (18)
< = i V е0,и |W;i<P" ' = 1-2--M. (19)
[w* (y - (A ) Xt), в противном случае , Фактически (A , z ) является оптимальным двойственным решением задачи (15)—(17).
Повышение эффективности алгоритма ВМНМ-оценивания
Известен точный алгоритм МНМ-оценивания [17]. Этот алгоритм имеет вычислительную сложность O(M2N2 + M4N log N + M2N log2 N) . Кусочно-линейная форма целевой функции и простая структура допустимого множества задачи (15)—(17) дает основание полагаться на возможность ее более эффективного решения задачи. Решения таких задач, как (15)—(17), возможны с помощью алгоритма, основанного на методе проекции градиента [8, 9]. В рассматриваемом здесь случае, благодаря простоте структуры допустимого множества, вычислительная сложность алгоритма не превышает величины O(M2N2) . Алгоритм PrGrad
Вход: XeRNxRM, ре( 1R+)M, jeMM
M
Выход: w* = arg max V wi ■ yi.
WE. , 1=1
Шаг 1. /* Инициализация*/
Шаг 1.1. k = 0 /* Счетчик итераций */
Шаг 1.2. w(0) = {w(0) = 0}I=12 м /* Координаты начальной точки */
Шаг 2. /* Текущая итерация к */
Шаг 2.1. Вычислить набор индексов ненасыщенных переменных:
S(к) = {i: - f < w(k) < pti = 1,2,., M}. Шаг 2.2. Вычислить подматрицы ненасыщенных переменных:
X(k) = X : i е S(k), j = 1,2,... N}, y(k) = {y; i е S(k)}. Шаг 2.3. Вычислить матрицу проектирования
p(k) = e - X(k)т ^X(k)X(k)T J-1 X(k)
и проекцию градиента целевой функции g(k) = P(k) y(k).
Шаг 2.4. Проверить необходимое и достаточное условие максимума:
если g(k) = 0, то перейти на Шаг 3. Шаг 2.5. Вычислить максимально допустимую длину шага в направлении g(k -1:
a = arg max {a: -p < w(k) + ag(k) < p, i е S(k)}.
a ^ '
Шаг 2.6. Вычислить следующую точку
w(k+1) = w(k) +ag(k), i е S(k). Шаг 2.7. Положить k = k +1 и перейти на Шаг 2.1. Шаг 3. Завершение алгоритма: w* = w(k).
Конец описания алгоритма PrGrad
Теорема. Алгоритм PrGrad корректно решает проблему (15)—(17). Его вычислительная сложность не превышает величины O(M2N2).
Доказательство. Каждая k-я итерация алгоритма состоит из допустимого перемещения из текущей точки w(k) в следующую точку w(k+1) (см. Шаг 2.6) в направлении проекции g(k) градиента y целевой функции (15) на пересечение множества решений системы уравнений (16) и множества решения системы активных ограничений из (17), которые могут быть определены
множеством r(k) = |/ : |Wi(K^ = pt | насыщенных переменных.
Более того, если выполнено необходимое условие для экстремума, т. е. из g(k) = 0, то для всех i е R(k) выполняются условия w(k)y > 0 . Поэтому градиент y не может быть представлен как неотрицательная линейная комбинация градиентов (т. е. внутренних нормалей) активных ограничений. Поэтому в нашем случае в соответствии с теоремой Куна-Такера необходимое условие экстремума также является достаточным, и тело шага 2 будет выполнено не более чем N раз.
Вычислительная сложность шагов 2.1, 2.4-2.6 не превышает величины O(M).
Матрица X(k) имеет N строк и не более M столбцов. Поэтому вычислительная сложность шага 2.2 не будет превышать O (MN).
Опишем более подробно оценку вычислительной сложности шага 2.3. Вычислительная сложность умножения (N хM) -матрицы на (M х N) -матрицу, т. е. вычисление матрицы в
квадратных скобках, не превышает величины O(MN2). Сложность полученной обратной
(N х N) -матрицы не превышает значения O(N3) .
Вычислительная сложность умножения (M х N) -матрицы на (N х N) -матрицу не превышает O(MN2). Вычислительная сложность умножения (M х N) -матрицы на (N х M) -матрицу также не превышает значения O(MN2) .
Таким образом, вычислительная сложность тела цикла не превышает O(MN2 + N3), а вычислительная сложность алгоритма не превышает O(M2N2) , поскольку M > N .
Теорема доказана.
Заключение
Идентификация квазилинейного уравнения (1) авторегрессии порядка m, имеющего n экзогенных переменных, по M отсчетам траектории возможна применением ОМНМ на основе связи между ОМНМ- и ВМНМ-моделями, позволяющей решить проблему определения ОМНМ-оценок посредством итеративной процедуры ОМНМ-оценивания с ВМНМ-оценками. При этом задача ВМНМ-оценивания представляет задачу линейного программирования, содержащую N = n + m переменных и 2M ограничений. С помощью предложенного алгоритма PrGrad данная задача решается за время не превосходящее O(M2N2).
Литература
1. Райбман, Н.С. Что такое идентификация? / Н.С. Райбман. - М.: Наука, 1970. - 119 с.
2. Gabisch, G. Business Cycle Theory / G. Gabisch, H.W. Lorenz. - Springer-Verlag, 1989. -250 p.
3. Davies, S. Inter-Firm Diffusion of Process Innovations / S. Davies // European Economic Review. - 1979. - Vol. 12, Iss. 4. - P. 299-317.
4. Grabec, I. Synergetics of Measurements, Prediction and Control. Springer Series in Synergetics, Book 68 / I. Grabec, W. Sachse. - Berlin-New York: Springer Verlag, 1997. - 458 p.
5. Panyukov, A.V., Stable Identification of Linear Autoregressive Model with Exogenous Variables on the Basis of the Generalized Least Absolute Deviation Method / A.V.Panyukov, Y.A. Mezaal // Вестник ЮУрГУ. Серия «Математическое моделирование и програмирование». - 2018. - Т. 11, № 1.- С. 35-43.
Панюков А.В., Мезал Я.А.
Параметрическая идентификация квазилинейного разностного уравнения
6. Panyukov, A.V. Stable Estimation of Autoregressive Model Parameters with Exogenous Variables on the Basis of the Generalized Least Absolute Deviation Method / A.V. Panyukov, Y.A. Mezaal // IFAC-PapersOnLine. - 2018. - Vol. 51, Iss. 11. - P. 1666-1669.
7. Panyukov, A.V. Linkage Between Wlad and Glad and its Applications for Autoregressive Analysis / A.V. Panyukov, I.A. Tetin, Y.A. Mezaal // Proceedings of the 4th International Conference «Information Technologies for Intelligent Decision Making Support (ITIDS'2016)». - 2016. -С.224-227.
8. Мину, М. Математическое программирование. Теория и алгоритмы / М. Мину. -М.: Наука, 1990. - 485 с.
9. Rosen, J.B. The Gradient Projection Method for Nonlinear Programming, part 1: linear constraints / J.B Rosen // Journal S.I.A.M. - 1960. - Vol. 8. - P. 181-217.
10. Гурин, Л.С. О состоятельности оценок метода наименьших квадратов / Л.С. Гурин // Математическое обеспечение космических экспериментов: сб. науч. тр. - М.: Наука, 1978. -С. 69-81.
11. Мудров, В.И. Методы обработки измерений: квазиправдоподобные оценки / В.И. Мудров, В.Л. Кушко. - М.: URSS, 2014. - 302 c.
12. Shestakov, A.L. The Theory of Optimal Measurements / A.L. Shestakov, A.V. Keller, G.A. Sviridyuk // Journal of Computational and Engineering Mathematics. - 2014. - Vol. 1, № 1. -С. 3-16.
13. Huber, P. Robust Statistics / P.Huber, E M. Ronchetti. - New York, Wiley, 2009. - 354 p.
14. Bloomheld, P. Least Absolute Deviations - Theory, Applications, and Algorithms / P. Bloomheld, W. Steiger. - Boston-Basel-Stuttgart: Birkhauser, 1983. - 351 p.
15. Pan, J. Weighted Least Absolute Deviations Estimation for ARMA Models with Infinite Variance / J. Pan, H. Wang, Y. Qiwei // Econometric Theory. - 2007. - Vol. 23, Issue 5. - P. 852-879.
16. Panyukov, A. Stable Parametric Identification of Vibratory Diagnostics Objects / A. Panyukov, A. Tyrsin // Journal of Vibroengineering. - 2008. - Vol. 10, no. 2. - P. 142-146.
17. Тырсин, А.Н. Точные алгоритмы реализации метода наименьших модулей на основе спуска по узловым прямым / А.Н. Тырсин, A. Азарян // Вестник БГУ. Математика. Информатика. - 2017. - № 4. - С. 21-32.
Поступила в редакцию 20 сентября 2019 г.
Bulletin of the South Ural State University Series "Mathematics. Mechanics. Physics" _2019, vol. 11, no. 4, pp. 32-38
DOI: 10.14529/mmph190404 PARAMETRIC IDENTIFICATION OF QUASILINEAR DIFFERENCE EQUATION A.V. Panyukov, Ya.A. Mezal
South Ural State University, Chelyabinsk, Russian Federation E-mail: paniukovav@susu.ru
Identification of the quasilinear difference equation is reduced to the problem of regression analysis with mutually dependent observable variables. This makes the classical solution schemes, based on the least squares method and its variations, ineffective. Finding estimates of the autoregressive equation coefficients is significantly complicated by poor conditionality of the system of equations, which represent necessary conditions for the minimum sum of squared deviations. In this case, estimates of the problem parameters are untenable. For solving such problems, it is possible to use the generalized least absolute deviations (GLAD) method, reduced to problems sequence of estimates of the autoregressive equation coefficients with the weighted least absolute deviations (WLAD) method. The article proposes an algorithm for solving the problem of WLAD-estimation, based on its conversion to the problem of linear programming (LP) of simple structure. The simplicity of the admissible set of the LP problem structure lies in the intersection of a linear subspace with a parallelepiped. It allows to propose an effective algorithm for its solution based on the gradient projection method. The algebraic computational
complexity of the proposed algorithm does not exceed the value O(N2M2), where N is the number of coefficients in the equation under study, and M is the number of the observed values. This WLAD computational complexity estimate is considered the best among the known ones.
Keywords: least absolute deviations; autoregressive model; linear programming; gradient projection method; computational complexity.
References
1. Raybman N.S. Chto takoe identifikatsiya? (What is Identification?). Moskow, Nauka Publ., 1970, 119 p. (in Russ.).
2. Gabisch G., Lorenz H.W. Business Cycle Theory. Springer-Verlag, 1989, 250 p. DOI: 10.1007/978-3-642-74715-1
3. Davies S. Inter-Firm Diffusion of Process Innovations. European Economic Review, 1979, Vol. 12, Iss. 4, pp. 299-317. DOI: 10.1016/0014-2921(79)90023-0
4. Grabec I., Sachse W. Synergetics of Measurements, Prediction and Control. Springer Series in Synergetics, Book 68, Berlin-New York, Springer Verlag, 1997, 458 p.
5. Panyukov A.V., Mezaal Y.A. Stable Identification of Linear Autoregressive Model with Exogenous Varia-Bles on the Basis of the Generalized least Absolute Deviation Method. Bulletin of the South Ural State University. Series "Mathematical Modelling, Programming & Computer Software ", 2018, Vol. 11, no. 1, pp. 35-43. DOI: 10.14529/mmp180104
6. Panyukov A.V., Mezaal Y.A. Stable estimation of autoregressive model parameters with exogenous variables on the basis of the generalized least absolute deviation method. IFAC-PapersOnLine, 2018, Vol. 51, Iss. 11, pp. 1666-1669. DOI: 10.1016/j.ifacol.2018.08.217
7. Panyukov A.V., Tetin I.A., Mezaal Y.A. Linkage Between Wlad and Glad and its Applications for Autoregressive Analysis. Proceedings of the 4th International Conference "Information Technologies for Intelligent Decision Making Support (ITIDS'2016)", 2016, C. 224-227.
8. Minoux M. Programmation mathematique: theorie et algorithmes. Paris, Dunod, 1983. (in French).
9. Rosen J.B. The Gradient Projection Method for Nonlinear Programming, Part 1: Linear Constraints. Journal S.I.A.M., 1960, Vol. 8, P. 181-217.
10. Gurin L.S. O sostoyatel'nosti otsenok metoda naimen'shikh kvadratov (On the Consistency of the Least Squares Method Estimates). Matematicheskoe obespechenie kosmicheskikh eksperimentov: sb. nauch. tr. (Mathematical Support of Space Experiments: Collection of Scientific Papers). Moscow, Nauka Publ., 1978, pp. 69-81. (in Russ.).
11. Mudrov V.I., Kushko V.L. Metody obrabotki izmereniy: kvazipravdopodobnye otsenki (Methods for Processing Measurements: Quasilike Estimates). Moscow, URSS Publ., 2014, 302 p. (in Russ.).
12. Shestakov A.L., Keller A.V., Sviridyuk G.A. The Theory of Optimal Measurements. Journal of Computational and Engineering Mathematics, 2014, Vol. 1, no. 1. pp. 3-16. (in Russ.).
13. Huber P., Ronchetti EM. Robust Statistics. New York, Wiley, 2009, 354 p. DOI: 10.1002/9780470434697
14. Bloomheld P., Steiger W. Least Absolute Deviations - Theory, Applications, and Algorithms. Boston-Basel-Stuttgart, Birkhauser, 1983, 351 p. DOI: 10.1007/978-1-4684-8574-5
15. Pan J., Wang H., Qiwei Y. Weighted Least Absolute Deviations Estimation for ARMA Models with Infinite Variance. Econometric Theory, 2007, Vol. 23, Issue 3, pp. 852-879. DOI: 10.1017/S0266466607070363
16. Panyukov A., Tyrsin A. Stable Parametric Identification of Vibratory Diagnostics Objects. Journal of Vibroengineering, 2008, Vol. 10, no. 2, pp. 142-146.
17. Tvrsin A.N., Azaryan A.A. Exact Algorithms for Implementation of the Least Absolute Deviations' Method Based on the Descent Through the Nodal Straight Lines. BSU Bulletin. Mathematics, Informatics, 2017, Issue 4, pp. 21-32. (in Russ.).
Received September 20, 2019
