Проекционный алгоритм линейного программирования на основе поиска пересечения прямой и зонотопа

Abstract

Аннотация. В работе рассматривается проекционный алгоритм решения задачи линейного программирования при наличии интервальных ограничений на все переменные. Ограничения задачи преобразуются в многогранник специального вида - зонотоп, а ее решение соответствует минимальной по некоторой координате точке пересечения зонотопа с прямой. Предложена параллелизация алгоритма на основе метода типа бисекции. Ключевые слова: линейное программирование, параллельные вычисления, выпуклые многогранники, зонотопы, проекция.

-

Copyright © by the paper's authors. Copying permitted for private and academic purposes.

In: A. Kononov et al. (eds.): DOOR 2016, Vladivostok, Russia, published at http://ceur-ws.org В настоящей работе задача ЛП рассматривается в контексте т.н. встроенной оптимизации (embedded optimization) – направления, развивающегося с начала этого столетия, по-видимому, исключительно за рубежом (см. напр. недавние англоязычные работы по автоматическому управлению только за один 2014 г. [ 13-17 ]) – автору неизвестны отечественные работы на эту тему. Встроенная оптимизация предполагает решение оптимизационных задач небольшого размера, но в режиме реального времени на базе маломощного (по сравнению с персональным, или, тем более, кластерным компьютером) встроенного микропроцессора для задач управления или обработки сигналов. В работах по встроенной оптимизации принимают участие такие известные зарубежные специалисты, как Стивен Бойд [ 18 ], Ю.Е. Нестеров, созданы системы генерации оптимизационных кодов для встроенных процессоров, такие как ACADO [ 19 ], FORCES [ 20 ], Multi-Parametric Toolbox [ 21 ].

При решении задач встроенной оптимизации возникают технические ограничения, которые могут показаться необычными исследователям, занимающимся задачами оптимизации применительно к персональным или тем более кластерным компьютерам. Так, например, во многих работах по встроенной оптимизации считается преимуществом полное отсутствие необходимости решения систем линейных уравнений, обращения матриц. Причиной является как отсутствие необходимых стандартных библиотек функций для целевого процессора с операционной системой реального времени или без нее, так и повышенного за счет этих операций времени выполнения программы. Гораздо более важной, чем в традиционных задачах математической оптимизации, становится незавышенная оценка сверху времени выполнения программы для всех исходных данных из некоего заданного диапазона, так как по окончании вычислений происходит выдача команд управления или обработка информации для технического объекта, которая регулярно повторяется с заданным периодом в диапазоне микро- и миллисекунд. Немаловажной является и простота реализации алгоритма с точки зрения инженера-специалиста в конкретной прикладной области, зачастую не являющегося экспертом в теории оптимизации.

Отметим, что в задачах управления реального времени весьма популярен вариант метода «быстрого градиента» Ю.Е. Нестерова [ 22 ], основанный на проекции на гиперкуб (которая выполняется покомпонентно путем выбора минимального из двух чисел). В работах [ 23,24 ] метод используется для решения задачи минимизации квадратичной положительно определенной функции с интервальными ограничениями на переменные, получены оценки сверху времени выполнения метода для заданной точности решения. В работе [ 25 ] метод используется в качестве вспомогательного для решения задачи квадратичного программирования с линейными ограничениями общего вида. Предложенные в этих работах подходы, безусловно, относятся к методам проекционного типа и представляют интерес и за пределами задач встроенной оптимизации. Так, например, отсутствие в алгоритме поиска решения системы линейных уравнений может оказаться выигрышным и для задач «гигабайтной» оптимизации [ 5 ] – как известно, квадратичное программирование широко используется в задачах машинного обучения, в частности в методе SVM (support vector machines) [ 26 ],

Мотивацией для настоящей работы послужили исследования в области использования полиэдральных (т.е. основанных на минимизации l1- и l∞-норм) критериев качества в задачах управления на основе предсказательного моделирования (model predictive control) [ 27-29 ] и другие подходы к управлению на основе ЛП в реальном времени (см. напр. [ 30-32 ]), а также весьма популярная в последнее время задача восстановления сильно разреженных сигналов на основе l1-нормы (compressed sensing) [ 33-37 ]. Все эти задачи сводятся к задаче ЛП, а их выполнение в большинстве случаев подразумевается на встроенных микропроцессорах.

Исторически первой задачей управления, решенной с помощью ЛП, явилась задача оптимального по быстродействию управления для дискретных линейных систем c ограничениями на управление [ 38 ]. Интересно отметить, что первая попытка использовать конечномерную оптимизацию в задачах управления с обратной связью связана также с линейным программированием и является международно признанным приоритетом отечественных ученых в лице А.И. Пропоя [ 39 ].

В настоящей работе предпринята попытка с одной стороны, использовать многоядерность встроенного процессора (так, например, в современных архитектурах ARM и Intel Atom доступны по крайней мере 4-8 ядер) в процессе решения задачи ЛП, и с другой – построить алгоритм таким образом, чтобы можно было легко вычислить нужное число итераций для заданной точности решения. При этом в настоящей работе предприняты пока только лишь начальные шаги в выбранном направлении. Предложен параллельный алгоритм решения задачи, но пока не выполнены необходимые для оценки его работоспособности вычисления, не выполнено сравнение с другими возможными подходами к решению задач ЛП на встроенных микропроцессорах (например, ADMM [ 17,40 ]). 2 Линейное программирование как поиск пересечения прямой и зонотопа Зонотопы  выпуклые многогранники, являющиеся аффинными проекциями многомерного куба [41]:

Z  {z  Rn : z  z0  Hx, || x ||  1}, x  Rm , n  m,|| x ||  im1,a...x,n | xi |. (1) Здесь матрица

H  [h1 | h2 | ... | hm ] Rnm содержит столбцыгенераторы. В задаче линейного программирования (ЛП) с минимизируемой целевой функцией pT z и ограничением z  Z решение легко определяется в замкнутом виде: (2) если i, pT hi  0 . В случае i, pT hi  0 решение представляет собой выпуклую оболочку некоторых граничных точек Z . Эта формула может быть использована в методе условного градиента (название предложено Б.Т. Поляком [42], в англоязычной литературе метод больше известен как Frank-Wolfe algorithm по именам его изобретателей Маргарет Франк и Филипа Вульфа), где на каждом шаге метода ищется линейное приближение выпуклой функции и затем его минимум на Z . В последнее время метод условного градиента привлек внимание исследователей в области машинного обучения, где он дает преимущества при решении больших задач квадратичной оптимизации – см. напр. работы [ 43-45 ].

В [ 46 ] был предложен подход к решению задачи ЛП вида

min cT x, Ax  b, || x ||  1, с помощью выпуклой оптимизации на зонотопе. Задачу можно переформулировать с учетом данного выше определения как (3) (4) (5) b  A  min , l( )     Z, H  cT  , z0  0.

  Теперь мы ищем минимальную по  точку пересечения прямой l() и зонотопа. Заметим, что путем добавления дополнительных переменных и ограничения их некоторыми конечными интервалами в подобном виде можно, по сути, сформулировать любую практическую задачу ЛП. Для иллюстрации подхода на рис. 1 изображены зонотоп и прямая l( ) , построенные для простейшей задачи ЛП: max x , x  x2  0, x1  x2  1, || x ||  2.

1 1 Для преобразования неравенства задачи в равенство вводится дополнительная переменная x3 [ 0,5 ] (здесь x3 =5 соответствует x1,2  2 ). Интервал [-2,2] для значений   x1 определяется ограничениями на x1 . Для преобразования в форму (1) столбцы H в (4) необходимо поделить на ширину интервалов. Тогда

0  1 / 4 l( )  1  , H   1 / 4   1 / 4 1 / 4 (6) В данном случае (как и в формуле (5)) зонотоп несимметричен (смещен) относительно начала координат, это соответствует ненулевому z0 в (4). Рис. 1. Зонотоп (зеленый) и прямая l() (красная), построенные для примера. Минимальное значение вертикальной координаты =-0.5 в точке пересечения соответствует максимуму переменной x1 =0.5 в задаче ЛП.

Описанный подход позволяет решать задачу проекционными методами. Метод LP-Newton, предложенный в [ 46 ], решает задачу за конечное число шагов (весьма похожий метод предложен позднее также в [ 47 ] для другой задачи, а в [ 48 ] метод [ 46 ] был распространен на общий случай ЛП с проекцией на конус). Поясним суть метода. Вначале выполняется проекция произвольной точки, принадлежащей прямой, на зонотоп, для чего задача сводится к поиску точки многогранника, ближайшей к началу координат (перенесенного в точку на прямой). Эта задача в работе [ 46 ] решается методом из [ 49 ] (заметим, что в отечественной литературе по негладкой оптимизации для аналогичной задачи предложен известный метод МДМ [ 50 ]). Затем строится плоскость, касательная к поверхности уровня квадратичной функции евклидова расстояния в точке проекции на зонотопе и находится точка пересечения этой плоскости с прямой (автор работы [ 46 ] проводит аналогии с методом Ньютона для гладких функций, где также ищется пересечение с касательной плоскостью, отсюда название метода). Затем уже эта точка берется в качестве начальной и все операции повторяются.

Отметим, что в отечественной литературе конечношаговый метод решения задачи ЛП на другой основе был впервые предложен, по видимому, в [51] (развит далее в [ 52 ]), в последнее время конечношаговый метод ЛП был предложен также в работах под рук. акад. Ю.Г. Евтушенко [ 53,54 ].

В настоящей работе рассматривается возможность распараллеливания метода пересечения зонотопа и прямой на других принципах. Предлагается строить алгоритм на основе деления отрезка прямой l() на несколько интервалов (пропорционально количеству параллельных ветвей программы) и одновременного поиска проекций границ этих интервалов на зонотоп (одним из доступных методов – например, методом МДМ, методом Вульфа [ 49 ], методом условного или «быстрого» градиента). Необходимым условием для начала работы такого алгоритма является поиск начальной внутренней точки зонотопа, для чего теми же методами решается задача

b  A  min yT y, y      T  x, || x ||  1,  min     max .

  c  (7) Предлагаемый мультисекционный (имеется в виду аналог бисекции в случае нескольких интервалов деления отрезка) алгоритм состоит из следующих шагов:

n n  Шаг 1. Вычислить  min   cisign(ci ),  max   cisign(ci ), i1 i1 что соответствует оптимизации сT x на ограичениях || x ||  1.  Шаг 2. Вычислить по (7) внутреннюю точку l( * ) и присвоить  max   * .  Шаг 3.

Исходя из числа N паралельных процессов, разбить интервал [ min ,  max ] на N равных подинтервалов. Для их граничных точек одновременно вычислить их проекции на зонотоп. Пометить точки, проекция которых не равна самим этим точкам, как внешние.  Шаг 4. Присвоить границам интервала новые значения – максимальное среди внешних точек и минимальное среди внутренних, если разница между соответствующими точками мала – выйти, иначе перейти на Шаг 3. Целью работы является предложить для дальнейшего исследования и обсуждения новый подход к решению задач ЛП на основе на основе параллельного мультисекционного деления отрезка и проекций на многогранник специального вида. Автор не утверждает, что рассматриваемый в работе метод ЛП имеет какие-либо преимущества перед предлагавшимися ранее в контексте задач ЛП большой размерности, в то же время, поскольку в методе отсутствует как таковая операция решения системы уравнений, представляют интерес дальнейшие исследования также и в этом направлении.

Заметим, что в работе описан лишь один из множества возможных методов, использующих общую идею поиска пересечения прямой и зонотопа. Очевидно, что решение может быть получено и классическим методом альтернирующих проекций [ 7,55,56 ], в котором поочередно находятся проекции текущей точки прямой l() на зонотоп и полученной точки на зонотопе обратно на прямую. В контексте рассматриваемых в работе подходов представляет интерес также работа [ 57 ], где проекционными методами решается задача смешанного целочисленного программирования без использования стандартного в таких случаях метода ветвей и границ.

Автор благодарен анонимному рецензенту, обратившему его внимание на работы отечественных авторов по конечношаговым методам ЛП, а также участникам конференций SIAM Conference on Parallel Processing for Scientific Computing (Париж, апрель 2016 г.) и International Congress on Mathematical Software (Берлин, июль 2016 г.), где проходило обсуждение вариантов метода. Список литературы

Maxim Demenkov Abstract. We study linear programming with box constraints on all variables (in addition to other linear constraints). Following research of Fujishige et al., we consider this problem as finding an intersection between a line and a specially constructed zonotope. This approach allows us to explore a variety of projection-type optimization algorithms, where we do not need to rely on the solution of linear equations.

References

1. Бердникова Е. А. , Ерeмин И . И., Попов Л . Д.: Распределенные фейеровские процес- сы для систем линейных неравенств и задач линейного программирования . Автома- тика и телемеханика, № 2 , с. 16 - 32 ( 2004 )

5. Нурминский Е.А. , Шамрай Н .Б.: Полиэдры, независимые множества и гигабайтное линейное программирование . Материалы V Всерос. конф. « Проблемы оптимизации и экономические приложения», Омск: Изд-во Омск. гос. ун-та, c . 52 - 56 ( 2012 )

6. Nurminski , E.A. : Single-projection Procedure for Linear Optimization . J. of Global Optimization, Online First Articles ( 2015 )

8. Ерeмин И. И. , Попов Л . Д.: Фейеровские процессы в теории и практике: обзор по- следних результатов . Изв. вузов. Матем ., № 1 , с. 44 - 65 ( 2009 )

9. Нурминский Е .А.: Фейеровские алгоритмы с адаптивным шагом . Ж. вычисл. матем. и матем . физ., № 5 , c. 791 - 801 ( 2011 )

12. Necoara I. , Nesterov Y. , Glineur F. : A Random Coordinate Descent Method on Large Optimization Problems with Linear Constraints . Technical Report , University Politehnica Bucharest ( 2011 )

13. Guiggiani A. , Patrinos P. , Bemporad A. : Fixed-point Implementation of a Proximal Newton Method for Embedded Model Predictive Control . In Proc. of the 19th IFAC World Congress, Cape Town , South Africa ( 2014 )

14. Rubagotti M. , Patrinos P. , Bemporad A. : Stabilizing Linear Model Predictive Control under Inexact Numerical Optimization . IEEE Trans. on Automatic Control 59 , 1660 - 1666 ( 2014 )

15. Korda M. , Jones C. N. : Certification of Fixed Computation Time First-order Optimizationbased Controllers for a Class of Nonlinear Dynamical Systems . In Proc. of the American Control Conference , Portland, Oregon ( 2014 )

16. Zeilinger M. N. , Raimondo D. M. , Domahidi A. , Morari M. , Jones C.N. : On Real-time Robust Model Predictive Control . Automatica 50 , 683 - 694 ( 2014 )

17. Jerez J.L. , Goulart P.J. , Richter S. , Constantinides G.A. , Kerrigan E.C. , Morari M. : Embedded Online Optimization for Model Predictive Control at Megahertz Rates . IEEE Trans. on Automatic Control 59 , 3238 - 3251 ( 2014 )

18. Wang Y. , Boyd S. : Fast Model Predictive Control using Online Optimization . IEEE Trans. on Control Systems Technology 18 , 267 - 278 ( 2010 )

19. Houska B. , Ferreau H. J. , Diehl M.: An Auto-generated Real-time Iteration Algorithm for Nonlinear MPC in the Microsecond Range . Automatica 47 , 2279 - 2285 ( 2011 )

20. Domahidi A. , Zgraggen A. , Zeilinger M.N. , Morari M. , Jones C.N. : Efficient Interior Point Methods for Multistage Problems Arising in Receding Horizon Control . In: Proc. of the CDC , 668 - 674 , Maui, USA ( 2012 )

21. Herceg M. , Kvasnica M. , Jones C.N. , and Morari M. : Multi-Parametric Toolbox 3.0 . In: Proc. of the European Control Conference , 502 - 510 , Zurich, Switzerland ( 2013 ) [Элек- тронный ресурс]: URL: http://control.ee.ethz.ch/~mpt (Дата обращения: 10 . 07 .16).

22. Нестеров , Ю.Е.: Метод решения задачи выпуклого программирования со скоростью сходимости O(1/k 2 ) . Докл. АН СССР , т . 269 , № 3, с. 543 - 548 ( 1983 )

23. Richter S. , Jones C. N. , Morari M. : Real-time Input-constrained MPC using Fast Gradient Method . In: Proc. IEEE CDC ( 2009 )

24. Richter S. , Jones C. , Morari M. : Computational Complexity Certification for Real-time MPC with Input Constraints Based on the Fast Gradient Method . IEEE Trans. on Automatic Control 57 , 1391 - 1403 ( 2012 )

25. Patrinos P. , Bemporad A. : An Accelerated Dual Gradient-Projection Algorithm for Embedded Linear Model Predictive Control . IEEE Trans. on Automatic Control 59 , 18 - 33 ( 2013 )

26. Vapnik , V. : The Support Vector Method of Function Estimation . In: Suykens, J.A.K. , Vandewalle , J . (eds) Nonlinear Modeling : Advanced Black-Box Techniques . chap. 3 , 55 - 85 . Kluwer Academic Publishers, Boston ( 1998 )

27. Bemporad A. , Borrelli F. , Morari M. : Model Predictive Control based on Linear Programming - the Explicit Solution . IEEE Trans. on Automatic Control 47 , 1974 - 1985 ( 2002 )

28. Borrelli F. , Bemporad A. , Morari M. : Predictive Control for Linear and Hybrid Systems . [Электронный ресурс]: URL: http://www.mpc.berkeley.edu/mpc-course-material (Дата обращения: 10.07 .16)

29. Филимонов Н .Б.: Методы полиэдрального программирования в дискретных задачах управления и наблюдения. Методы классической и современной теории автоматиче- ского управления. Учебник в 5-и тт . Т . 5. Методы современной теории автоматиче- ского управления . Гл. 7 . М.: Изд-во МГТУ им . Н.Э . Баумана, c. 647 - 720 ( 2004 )

30. Lazar M. , Jokic A. : Synthesis of Trajectory-dependent Control Lyapunov Functions by a Single Linear Program . In: Tabuada P., Majumdar R . (Eds.). Proceedings HSCC 2009 , 237 - 251 , LNCS No. 5469 . Berlin: Springer ( 2009 )

31. Нгуен Х.Н. , Гутман П .О., Олару С. , Ховд М .: Управление с ограничениями для ли- нейных стационарных систем: интерполяционный подход . Автоматика и телемеха- ника, № 1 , c. 68 - 89 ( 2014 )

32. Filimonov N.B. : Barrier Regulator Design using Polyhedral Predictive Control . In: Proc. of the IEEE International Conference “Stability and Control Processes” in memory of V.I. Zubov (SCP) , 48 - 51 , St. Petersburg ( 2015 )

33. Акимов П.А. , Деревянкин А .В., Матасов А .И.: Гарантирующий подход и l1- аппроксимация в задачах оценивания параметров БИНС при стендовых испытаниях . М.: Изд-во Московского университета ( 2012 )

34. Bartels R.H. , Conn A.R. , Sinclair J.W. : Minimization Techniques for Piecewise Differentiable Functions: the l1-solution to an Overdetermined Linear System . SIAM J. Numer. Anal . 15 , 224 - 241 ( 1978 )

35. Donoho D.L. For Most Large Underdetermined Systems of Equations, the Minimal l1- norm Solution is Also the Sparsest Solution . Communications on pure and applied mathematics , 59 , 797 - 829 ( 2006 )

36. Pudlewski S. , Melodia T. , Prasanna A. , Compressed-sensing-enabled Video Streaming for Wireless Multimedia Sensor Networks . IEEE Trans. Mobile Computing 11 , 1060 - 1062 ( 2012 )

37. Li S. , Xu L.D. , Wang X . Compressed Sensing Signal and Data Acquisition in Wireless Sensor Networks and Internet of Things . IEEE Trans. on Industrial Informatics 9 , 2177 - 2185 ( 2013 )

38. Zadeh L. , Whalen B. : On Optimal Control and Linear Programming . IRE Trans. on Automatic Control 7 ( 4 ), 45 - 46 ( 1962 )

39. Пропой А .И.: Применение методов линейного программирования для синтеза им- пульсных автоматических систем . Автоматика и телемеханика, № 7 , c. 912 - 920 ( 1963 )

40. Boley D. : Local Linear Convergence of the Alternating Direction Method of Multipliers on Quadratic or Linear Programs . SIAM J. Optim . 23 , 2183 - 2207 ( 2013 )

43. Beck , A. , Teboulle , M. : A Conditional Gradient Method with Linear Rate of Convergence for Solving Convex Linear Systems . Math. Methods of Op. Res . 59 , 235 - 247 ( 2004 )

44. Jaggi , M. : Revisiting Frank-Wolfe: Projection-free Sparse Convex Optimization . In: JMLR Proceedings 28 , 427 - 435 . ICML, Atlanta ( 2013 )

45. Lacoste-Julien , S. , Jaggi , M. : On the Global Linear Convergence of Frank-Wolfe Optimization Variants . In: Advances in Neural Information Processing Systems 28. NIPS , Montreal ( 2015 )

46. Fujishige , S. , Hayashi , T. , Yamashita , K. , Zimmermann , U. : Zonotopes and the LPNewton Method . Optimization and Engineering 10 , 193 - 205 ( 2009 )

47. Helmling , M. , Ruzika , S. : Towards Combinatorial LP Turbo Decoding . In: IEEE International Symposium on Information Theory , pp. 1491 - 1495 . IEEE Press, Istanbul ( 2013 )

48. Kitahara , T. , Mizunoa , S. , Shi , J.: The LP-Newton method for standard form linear programming problems . Operations Research Letters 41 , 426 - 429 ( 2013 )

49. Wolfe P. : Finding the Nearest Point in a Polytope . Mathematical Programming 11 , 128 - 149 ( 1976 )

50. Митчелл Б. Ф. , Демьянов В . Ф., Малоземов В . Н.: Нахождение ближайшей к началу координат точки многогранника . Вестник ЛГУ ., № 19 , c. 38 - 45 ( 1971 )

52. Лебедев В .Ю.: Декомпозиционный метод решения блочных задач линейного про- граммирования со связывающими переменными . Ж. вычисл. матем. и матем . физ.,, № 4, c. 881 - 886 ( 1981 )

54. Evtushenko Y. , Golikov A. : Linear Programming Projection Algorithms . Wiley Encyclopedia of Operations Research and Management Science ( 2010 )

55. Bauschke , H.H. , Borwein , J.M.: Dykstra's Alternating Projection Algorithm for Two Sets . J. Approx. Theory 79 , 418 - 443 ( 1994 )

56. Escalante , R. , Raydan , M. : Alternating Projection Methods . SIAM, Philadelphia ( 2011 )

57. Nowak I. : Column Generation-based Alternating Direction Methods for Solving MINLPs . [Электронный ресурс]: URL: http://www.optimization-online.org/DB_HTML/ 2015 /12/ 5233.pdf (Дата обращения: 10 . 07 .16)