CC BY
54
УДК 519.853
РЕШЕНИЕ ЗАДАЧИ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ С ИСПОЛЬЗОВАНИЕМ ОПЕРАТОРА-ПРОЕКТОРА
О.Н. Вылегжанин, Г.И. Шкатова
Институт «Кибернетический центр» ТПУ E-mail: [email protected]
Предложен метод решения задачи линейного программирования, основанный на вычислении оператора-проектора на пространство векторов активных ограничений. Оператор-проектор вычисляется посредством процедуры рекуррентного псевдообращения, что обеспечивает более высокую устойчивость вычислений по сравнению с преобразованием Гаусса-Жордана, используемого в симплекс-методе. Метод позволяет в рамках единой процедуры учесть наличие ограничений-равенств, вырожденность матрицы ограничений неравенств.
Ключевые слова:
Линейное программирование, оператор-проектор, рекуррентное псевдообращение.
Для решения задачи линейного программирования традиционно используют симплекс-метод Данцига, в основе которого лежит одна из форм преобразования Гаусса-Жордана [1]. Предлагаемая нами процедура решения задачи обеспечивает более высокую устойчивость вычислений по сравнению с преобразованием Гаусса-Жордана и позволяет в рамках единой процедуры учесть наличие ограничений-равенств и вырожденность матрицы ограничений неравенств, повышает точность получаемого решения и, кроме того, существенно сокращает объем вычислений.
Постановка задачи
Найти минимум функции: f =< g,x >| A-x< b,
(1)
максимального объема [3], получим матрицу полного ранга. Обозначим т - ранг матрицы А.
Для обоснования алгоритма получения опорного решения докажем два утверждения:
Утверждение 1. Если х - начальная точка, а ближайшая к ней плоскость многогранника М имеет направляющим вектором а, то х - проекция начальной точки на эту плоскость определяется как:
Ь; — аТ ■ х0
xi — X— I
a t - a t
Доказательство. Преобразуем (3) к виду: b — aT - xq
t
a i - a ,
(3)
(4)
где § - заданный, а л - подлежащий определению векторы, А - матрица известных коэффициентов системы ограничений, Ь - вектор правых частей системы неравенств.
Стратегия поиска
Если в исходной задаче определены ограничения-равенства вида А1х=Ь1, то для их учета можно воспользоваться методом, описанным в [1].
Условия Куна-Таккера для задачи (1) имеют вид:
1. Ах* < Ь;
2. Яг (Ь — а ■ х*) = 0, 1= 1,...,п; (2)
3. g + АТЯ = 0,
где х* - решение задачи, Я - вектор неположительных коэффициентов (коэффициентов Лагранжа), аI - строки матрицы ограничений А. Как известно, эти условия необходимы и достаточны для существования экстремума в точке х* [2].
Множество значений х, удовлетворяющих ограничениям из (1), составляет многогранник М. Будем полагать, что А - матрица полного столбцового ранга. В противном случае, применив процедуру выбора из матрицы А базисного набора столбцов
Так как —————
T
а. - а.
I I
■ скаляр, то направление из
Хо на точку х1 совпадает с направляющим вектором плоскости а;, т. е. нормально к этой плоскости. Умножим обе части (3) слева на а?:
¿T
T T . — а, ' Х- T 1
а. • x —а. • x0 +——р-- - — • — — — ,
а, - а.
(5)
т. е. х1 принадлежит плоскости, заданной направляющим вектором а.
Утверждение 2. Если точка хк принадлежит пересечению к плоскостей, заданных направляющими векторами а^=1...к, то ее проекция на пересечение плоскости, заданной строкой а,, с пересечением плоскостей, заданных строками а|=1...к, составляющими матрицу Ак, есть:
Ь. — аТ ■ х к
хк+1 = хк --ап (6)
а ; 2 ■а;
где 0=(/-Ак-Ак+) - проектор на пространство, перпендикулярное пространству, натянутому на строки матрицы Ак, а Ак+ - матрица, псевдообратная к матрице Ак, I — единичная матрица.
Доказательство.
1. Преобразуем (6) к виду:
Ъ , — а,
Х1- , 1 х>
■О,- а.
• О • а,.
Так как Ь' а' Хк
а : • О-а,
- скаляр, то направление из
точки хк на хк+1 перпендикулярно подпространству, натянутому на строки матрицы Ак.
2. Умножим обе части (6) слева на а,т:
Ъ, — а7:
а7: • О • а,
•а • О • а = Ъ.
Следовательно, точка хк+1 принадлежит плоскости с направляющим вектором а.
3. Умножая (6) на матрицу Ак слева, получим: Ъ, — а: -
Ак • Хк+1 = Ак • Хк + "
• О • а,
-Ла = А -Хк
т. к. Оаг по определению принадлежит подпространству, перпендикулярному пересечению плоскостей, задаваемых строками матрицы Ак, то Ак-О.а;=0.
Утверждения 1 и 2 обосновывают процесс получения опорного решения задачи линейного программирования, т. е. получения точки, принадлежащей вершине многогранника М. Выбрав произвольную начальную точку хо, выполнив (3), перейдем к следующей точке на плоскости, ближайшей к хо, а затем выполним последовательно (6) для тех г, для
Ъ — а: • х
которых у ==Ь-а;х,<0, а —к- - минимален. В
а : •а1
результате получим точку, для которой у<0| г=1...и, т. е. являющуюся вершиной многогранника М.
Получение оптимального решения связано с перемещением рабочей точки из этой вершины в ту, для которой значение целевой функции оптимально. Для шага перемещения необходимо из набора к ограничений-неравенств (строк матрицы Ак), которые обращаются в равенства в этой вершине, исключить одно неравенство и добавить другое. Соответственно, необходимо пересчитать новую матрицу Ак+. Для пересчета можно использовать процедуру рекуррентного псевдообращения [4]. Обозначим Ск матрицу, содержащую к столбцов, а ск - добавляемый столбец.
Тогда [4]:
ск = (Ск—1. ск )>
С =
к
С.
к—1[1 — ск• К
где к, =
С1 <1— Ск—1'С1) С\ •(1 — Ск—1 •Ск+—Ск
А также [5]:
Ск—1 = Нк—г( 1—ак к),
Вычислим коэффициенты Лагранжа Х=(Аку^. Обозначим / - множество индексов соответствующих неравенствам, включенным в Ак+. Для выбора неравенства а, исключаемого из / на первом шаге, нужно взять максимальное среди положительных значений коэффициентов Лагранжа.
Номер неравенства, включение которого в / соответствует новой вершине, соответствует минимальному положительному значению коэффициентов Лагранжа:
Ь , — а, • хк
а =--,
а , • О^ а,
Решению задачи соответствует точка х*, для которой все коэффициенты Лагранжа Я<0.
Продемонстрируем предложенный метод на следующем численном примере.
Пример. Требуется найти минимум функции у = / (х) = х1 — 2 • х2 + х3 — 3,424 • х4 + 0,5 • х5 — 6,25
при ограничениях:
А-х=Ь и Ах<Ь,
где:
А =
1,2 1,4 —0,7 0,26
—3 —1,5 0,5
А =
0, 4 —2
—1 0,1
—0,5 0,1
4,5 —3
4,5 —3
—1 2
2, 4 —1
—1 1
—0,45 0,2 0,1 —2 0,5 —10,1 0
0,5 —3,5 —3
0 —1 —1 —0,1 0,501 9,994 —7,481 —16,551 —2,398 —1,2 4
4,687
—1 —1 1
—0,2 0,35 —3,575 5, 425 —30,05
—5 —0,35 —5,2 —1,25
—0,5
0,5
, Ь =
4,5
0
2, 225
11,037
6,537
110,225
, Ь =
—12
—1,125
0,2
—2,625
Решение задачи поиска экстремума разбивается на три подзадачи:
1. Учет ограничений равенств.
2. Определение ранга матрицы, полученной преобразованием матрицы коэффициентов.
3. Поиск экстремума преобразованной задачи.
Для учета ограничений-равенств выполним следующие шаги:
1. Выберем столбцы матрицы А, образующие базисный набор. Это столбцы: 2, 4, 5. Ранг матрицы А равен трем: (а1т-Л3-а1)<1-10-5.
2. Найдем матрицу, псевдообратную к выделенным столбцам:
497549 0,252451 —0,250000 —0,024510 А
А+ =
где йк=ккт/кк\т.
0,253676 0,502451
—0,621324 —0,747549
—0,375000 0,750000
—0,036765 —0,024510
3. Умножим ее справа на ограничения-равенства:
А+ ■ А —
(1,00000 0,00000 0,00000 -0,07500 -0,50000 ^ 0,00000 1,00000 0,00000 -0,31250 -0,75000 0,00000 0,00000 1,00000 -0,97500 0,50000
Ч ' ' ' ' '
А+■ Ь — (-0,7500 -2,1250 2,7500)г.
4. Минимизируемая форма преобразуется к виду:
у — А + ■ / (X ) = = А+ ■ (х1 - 2 ■ х2 + х3 - 3,424 ■ х4 + 0> х5 - 6,25) = — 3 ■ Х1 Х3,
а ограничения-неравенства примут вид А'-х<Ь', где
А —
-0,998 8 -1 -3 2 16
-7 5 4,999
-26,998 2,998 -25 -5 , ь — 135 15
-0,162 0,4 1
0,455 -3 9,5
2 0,5 7
Матрица А' - полного столбцового ранга
(гкА=2).
Для поиска экстремума преобразованной задачи (4, 5) выполняем следующие шаги:
1. Нумеруем неравенства, определенные в (5), в порядке 0, 1, 2, ..., 7.
2. Берем начальную точку х0=(2; -8)т.
3. Вычисляем у=Ь—А'Х0.
4. На основании значений у из неравенств, определенных в (5), находим неравенство, ближайшее к точке х0. Это неравенство с номером 3.
5. Вычисляем проекцию на линию, определенную неравенством с номером 3. Это точка х1=(2,219; -7,797)т.
6. Вычисляем для новой точки у=Ь—А'Х{:
у — (7,578 25,145 -59,519 0 47,326 -4,474 14,912 -6,4 60)Г.
7. Точка х лежит на прямой, определенной неравенством с номером 3, а минимальное положительное значение соответствует неравенству с номером 0.
8. Вычисляем точку пересечения линий, определенных неравенствами с номерами 0 и 3: х2=(2; 0,0)т.
9. Вычисляем для новой точки у=Ь-А'Х2:
у — (0,000 0,000 9,000 -1,320 -18,999 -3,000)
Номера неравенств: 1, 0, 4, 5, 2, 7; исключено 3-е неравенство.
10. По вычисленным значениям определяем, что точка х3 принадлежит двум прямым, определенным неравенствами с номерами 0 и 1. Наличие положительных значений в у указывает, что опорное решение не достигнуто. Минимальное положительное значение (единственное) соответствует неравенству с номером 4.
11. Вычисляем точку пересечения линий, определенных неравенствами с номерами 4 и 1: х3=(2,5510204; 1,4693878)т.
12. Вычисляем для новой точки у=Ь—А'Х3
у — (0,000 0,000 -0,820 -15,509 -1,163)
Номера неравенств: 1, 4, 5, 2, 7; исключено 0-е неравенство.
13. Точка принадлежит прямым, определенным неравенствами с номерами 4 и 1. Среди элементов у нет положительных, т. е. опорное решение достигнуто. Проверим, является ли опорное решение оптимальным?
14. Вычисляем коэффициенты Лагранжа: Я=(-0,245; 0,347)т. Решение не оптимально, т. к. имеется положительное значение, соответствующее неравенству номер 4.
15. Исключаем неравенство с номером 4.
16. Определяем неравенство для включения - это неравенство номер 7.
17. Определяем точку пересечения 7 и 1: х4=(2,9; 2,4)т.
18. Вычисляем коэффициенты Лагранжа: А=(0,05; -1,7)т. Положительное значение соответствует неравенству номер 1.
19. Исключаем неравенство номер 1.
20. Вычисляем у=(-0,000 -13,299 -0,504)т.
21. Номера неравенств: 7, 2, 5.
22. Включаем неравенство 6. Точка пересечения: х5=(2,61; 3,545)т.
23. Вычисляем коэффициенты Лагранжа: А=(-1,545; -0,568)т. Решение оптимальное. Таким образом, предлагаемый метод позволяет
существенно сократить размерность решаемой задачи путем учета ограничений равенств, а также вырожденности матрицы ограничений-неравенств. Использование псевдообратной матрицы повышает точность вычисления коэффициентов Лагранжа и координат вычисляемой вершины многогранника. Применение операторов-проекторов гарантирует попадание из начальной точки в одну из вершин многогранника ограничений-неравенств, а затем удержание последовательности рабочих точек на поверхности этого многогранника. Приведенный численный пример подтверждает работоспособность метода.
СПИСОК ЛИТЕРАТУРЫ
1. Вылегжанин О.Н., Шкатова Г.И. Учет ограничений равенств при решении оптимизационных задач с линейными ограничениями // Известия Томского политехнического университета. - 2008. - Т. 312. - № 5. - С. 76-78.
2. Васильев Ф.П., Иваницкий А.Ю. Линейное программирование. - М.: Факториал, 1998. - 323 с.
3. Вылегжанин О.Н., Шкатова Г.И. Сравнительная оценка двух методов выбора наилучших линейных регрессоров // Применение математических методов и ЭВМ в медико-биологиче-
ских исследованиях / Межвузовский научно-техн. сб. под ред. В.А. Кочегурова. - Томск: Изд-во Томск. политехн. ин-та, 1988. - С. 18-22.
4. Гантмахер Ф.Р. Теория матриц. - 5-е изд. - М.: Физматлит, 2004. - 559 с.
5. Cline R.E. Representation for the Generalised Inverse of a Partitioned matrix // J. Soc. Industr. and Appl. Mathem. - 1964. - V. 12. -№ 3. - P. 588-600.
Поступила 15.04.2009 г.
УДК 681.31:533.95
МЕТОД СПЛАЙНОВОЙ ИНТЕРПОЛЯЦИИ ДЛЯ РЕШЕНИЯ УРАВНЕНИЯ ВЕКТОРНОГО ПОТЕНЦИАЛА В ЗАДАЧЕ ТРАНСПОРТИРОВКИ ЭЛЕКТРОННЫХ ПУЧКОВ
В.П. Григорьев, И.Л. Звигинцев, А.В. Козловских
Томский политехнический университет E-mail: [email protected]
Предложен метод сплайновой аппроксимации второй производной в узлах дискретной сетки по пространственной координате для решения уравнений векторного потенциала в электродинамических задачах при переходе от систем уравнений в частных производных к системам обыкновенных нелинейных дифференциальных уравнений. Показано, что разработанный метод сплайновой аппроксимации для такого рода задач работает на треть бы/стрее встроенных в Matlab стандартных функций сплайнов при относительно небольшом количестве узлов на сетке и равной погрешности.
Ключевые слова:
Сплайн-интерполяция, сплайновая аппроксимация, уравнения векторного потенциала, математическое моделирование.
В современных технологиях широко применяют низкоэнергетические сильноточные электронные пучки, как уникальные источники энергии. Для эффективного использования энергии, запасаемой в электронных пучках, требуется оптимизация их параметров и, следовательно, условий формирования и транспортировки на мишень. Для достижения этой цели из-за сложности физических процессов и дорогостоящих экспериментов большую роль приобретает численное моделирование и разработка новых математических методов [1, 2].
В данной работе предлагается новый метод для решения электродинамических задач. Апробация и тестирование этого метода проводится на примере решения проблемы транспортировки пучков к мишени в газе низкого давления.
Обычно подобного рода задачи решаются в пакете МаИаЬ, если другие пакеты с готовыми встроенными алгоритмами (например, Сош8о1) не могут справиться с поставленной задачей из-за расходимости вычислений. МаНаЬ имеет достаточно большой набор солверов, основанных на различных численных методах для решения систем обыкновенных дифференциальных уравнений. Для решения жестких систем с высокой точностью подходит солвер, основанный на многошаговом методе Гира, который допускает изменение порядка точности [3]. При вычислении производных в МаНаЬ ис-
пользуют сплайны. Для повышения скорости вычислений предложен метод, который работает быстрее, чем стандартные функции сплайн-интерполяции МаНаЬ на относительно небольшом количестве узлов сетки. При этом большего увеличения скорости вычислений можно достичь, решая данную задачу на кластере.
Исходная задача математически представлена системой нелинейных уравнений в частных производных:
dAz
и
ef
dnt ~dt
dt
1 A
r dr
d2A _ 4n
dr2
_— evbnb
(1)
_ G, ■ V. ■ И ■ П. + П ■ K ■ П
где А(г,1) - векторный потенциал поля, и;(г,0 -плотность ионов, ие(г,/) - плотность плазмы, г0 -классический радиус электрона, цДг,/) - частота соударений электронов и ионов плазмы, КС*,0 -коэффициент ионизации, - плотность газа, нь(г,1) - плотность пучка, уъ - скорость пучка, е -заряд электрона, ст; - сечение ионизации.
Эти нелинейные уравнения связаны через сложную зависимость приведенных выше параметров от Аг(г,/) и и;(г,0 [1].
Краевые условия:
