Решение задачи линейного программирования с использованием оператора-проектора Текст научной статьи по специальности «Математика»

УДК 519.853

РЕШЕНИЕ ЗАДАЧИ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ С ИСПОЛЬЗОВАНИЕМ ОПЕРАТОРА-ПРОЕКТОРА

О.Н. Вылегжанин, Г.И. Шкатова

Институт «Кибернетический центр» ТПУ E-mail: onv@am.tpu.ru

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

Ключевые слова:

Линейное программирование, оператор-проектор, рекуррентное псевдообращение.

Для решения задачи линейного программирования традиционно используют симплекс-метод Данцига, в основе которого лежит одна из форм преобразования Гаусса-Жордана [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: 8b20@sibmail.com

Предложен метод сплайновой аппроксимации второй производной в узлах дискретной сетки по пространственной координате для решения уравнений векторного потенциала в электродинамических задачах при переходе от систем уравнений в частных производных к системам обыкновенных нелинейных дифференциальных уравнений. Показано, что разработанный метод сплайновой аппроксимации для такого рода задач работает на треть бы/стрее встроенных в 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].

Краевые условия: