CC BY
224
REFERENCES
1. Gnedenko B.V., Kovalenko I.N. Vvedenie v teoriyu massovogo ob-sluzhivaniya [Introduction into queuing theory]. Moscow, Kom-Kniga Publ., 2007. 336 p.
2. Koenig D., Shtoyan D. Metody teorii massovogo obsluzhivaniya [Queuing theory methods]. Moscow, Radio i svyaz Publ., 1981. 128 p.
3. Goshtoni G. Sravnenie vychislennykh i modelirovannykh resulta-tov dlya puchkov soedinitelnykh liniy pri nalichii povtornykh po-pytok ustanovleniya svyazi [Comparison of calculated and simulated results for connecting line bunches at repeated attempts to establish communication]. VIIIITC Materials. Sydney, 1977, no. 1, pp. 1-16.
4. Falin G.I. Templeton J.G.C. Retrial queues. London, Chapman & Hall, 1997. 328 p.
5. Falin G.I. Asymptotic investigation of fully available switching systems with high repetition intensity of blocked calls. Moscow University Mathematics Bulletin, 1984, vol. 39, no. 6, pp. 72-77.
6. Falin G.I. A Survey of Retrial Queues. Queuing Systems,1990, vol. 7, pp. 127-167
7. Artalejo J.R., Gomez-Coral A. Retrial queuing systems: and computational approach. Berlin, Springer, 2008. 267 p.
8. Artalejo J.R., Joshua V.C., Krashnamoorthy A. An M/G/l retrial queue with orbital search by the server. Advances in Stochastic Modeling. New Jersey, Notable publications, 2002. pp. 41-54.
9. Nazarov A.A., Sudyko E.A. Metod asimptoticheskikh seminvari-antov dlya issledovaniya matematicheskoy modeli seti sluchayno-go dostupa [Method of asymptotic seven-invariants for research of mathematical model of a network of casual access]. Problemy peredachi informatsii, 2010, no. 1, pp. 94-111.
10. Nazarov A.A., Moiseyev S.P. Metody asimptoticheskogo analiza v teorii massovogo obsluzhivaniya [Methods of asymptotic analysis in the queuing theory]. Tomsk, NTL Publ. house, 2006. 112 p.
11. Borovkov A.A. Asimptoticheskie metody v teorii massovogo ob-sluzhivaniya [Asymptotic methods in queuing theory]. Moscow, Nauka Publ., 1980. 381 p.
УДК 519.688:622.276.5.001.42
АДАПТИВНАЯ ИНТЕРПРЕТАЦИЯ КРИВОЙ ВОССТАНОВЛЕНИЯ ДАВЛЕНИЯ ГОРИЗОНТАЛЬНЫХ СКВАЖИН C ДИАГНОСТИКОЙ ПОТОКОВ
Е.В. Романова, В.Л. Сергеев
Томский политехнический университет E-mail: Romanova88EV@mail.ru; SergeevVL@ignd.tpu.ru
Актуальность работы обусловлена необходимостью диагностики потоков при адаптивной интерпретации нестационарных гидродинамических исследований горизонтальных скважин с использованием аналитических моделей кривой восстановления давления.
Цель работы: разработка метода адаптивной интерпретации кривой восстановления давления, позволяющего диагностировать потоки в процессе проведения гидродинамических исследований горизонтальных скважин с одновременной оценкой фильтрационных параметров пласта.
Методы исследования: использованы теоретические и практические разработки в области гидродинамических исследований скважин, системного анализа, идентификации систем с учетом дополнительной априорной информации, оптимизации функций и линейной алгебры. Решение задач диагностики потоков, идентификации и интерпретации кривой восстановления давления проводилось на основе промысловых данных забойного давления на скважине с учетом экспертных оценок фильтрационных параметров пласта с использованием компьютерной программы Saphir.
Результаты: разработан адаптивный метод интерпретации кривой восстановления давления с диагностикой радиального и линейного потоков с одновременной оценкой фильтрационных параметров нефтяного пласта в процессе гидродинамических исследований горизонтальных скважин. На примере обработки результатов гидродинамических исследований по кривой восстановления давления двух горизонтальных скважин однородно-пористого нефтяного пласта показано, что метод адаптивной интерпретации позволяет повысить надежность определения времени начала радиального и линейного потоков.
Ключевые слова:
Интерпретация, диагностика, гидродинамические исследования, анализ нефтяных скважин, априорная информация, горизонтальные скважины, нефтяные пласты.
Введение
Известно, что интерпретация кривой восстановления давления (КВД) вертикальных скважин существенно отличается от интерпретации горизонтальных скважин, где сложный пространственный поток частиц жидкости сведен в определенные моменты времени к плоским фильтрационным потокам, представленным соответствующими уравнениями [1—3]. Так, например, радиальный поток горизонтальной скважины представлен уравнением забойного давления Рз(^ вида
CquB (2,25k t^
AP3(t) = P(t) - P(to) = ~r^~lnI----f- I, (1)
I m Mrnp )
т цгщ
а время его начала tr определяется по формуле
Гг = 1800^ 2т^С / к2, (2)
где Pз(t0) - забойное давление в момент остановки скважины; ^=л1 кк - радиальная проницаемость; кг, ку - вертикальные и горизонтальные проницаемости; q - дебит скважины перед ее остановкой, [№ - вязкость и объемный коэффициент нефти со-
ответственно; т - пористость; гпр - приведенным радиус скважины; Ь - длина горизонтальной части ствола скважины; d - расстояние от горизонтального ствола скважины до ближайшей границы нефтяного пласта; С - общая сжимаемость системы нефти и скелета породы; Се - константа, зависящая от используемой системы единиц, например, в системе СИ С8=9,205.
Линейный поток горизонтальной скважины представлен уравнением
АР =
шку
4пк1уЬ
Б
с временем его начала
ґ, = 160 І2ш^С / ку
(3)
(4)
На рис. 1 для скважины № 2031 достаточно четко виден прямолинейный участок, представляющий радиальный поток, начало и конец горизонтального участка, определяющий линейный поток с наклоном производной забойного давления порядка 45 градусов, что позволяет корректно применить известные методы интерпретации с использованием моделей забойного давления (1), (3) [1-3]. На рис. 2 для скважины № 1434 наблюдается иная ситуация, показывающая, что четко выделить радиальный и линейный потоков не представляется возможным, и возникает необходимость использования количественных критериев.
где ку - горизонтальная проницаемость; 5 -скин-фактор скважины при линейной геометрии потока.
Следует отметить, что на практике при исследовании скважин по КВД используются преимущественно эти два потока для сокращения времени простоя скважин.
Особенность интерпретации горизонтальных скважин с использованием аналитических моделей КВД заключается в выделении на КВД типа потока и решении обратной задачи идентификации по определению фильтрационных параметров пласта и скважины с использованием моделей (1), (3).
Рассмотрим два часто используемых способа выделения потоков. Первый, аналитический способ, заключается в определении времени начала потока по приведенным в (2) и (4) формулам. Второй, часто используемый на практике графо-ана-литический способ, основан на качественном критерии анализа производной от забойного давления Рз^)=йРз(^/Л в координатах ^Рз^)-^^) [1, 4, 5].
Для пояснения проблемы выделения потоков с использованием производной от забойного давления на рис. 1, 2 приведены значения забойных давлений и их производные в координатах ^Рз^)-^(^ для горизонтальных скважин однородно-пористого пласта нефтяного месторождения Тюменской области.
0,1 1 Время (ч)
Рис. 2. Кривая восстановления давления и ее производная горизонтальной скважины № 1434
Однако проблемным моментом при использовании аналитической оценки времени начала потоков является присутствие в соответствующих формулах (2), (4) фильтрационных параметров пласта, вертикальной и горизонтальной проницаемости к, ку, которые нам изначально неизвестны и подлежат определению.
Для решения указанных выше задач определения времени начала потоков в работе предлагается и исследуется метод адаптивной интерпретации КВД, позволяющий диагностировать потоки в процессе проведения гидродинамических исследований горизонтальных скважин с одновременной оценкой фильтрационных параметров пласта.
Диагностика потоков
по методу адаптивной интерпретации
Основу метода диагностики потоков представляет критерий вида
Д/* = а^шт3(Р* - р(«Ж)),
Аґ = [ґ + т] є[ґо, ґп]),
(5)
0,1 1 Время (ч)
Рис. 1. Кривая восстановления давления и ее производная горизонтальной скважиныы № 2031
где запись а^ тіпДх) означает точку минимума х* функции /(х) (Дх*)=тіпДх)); J - показатель качества заданный в видхе известной функции (либо функционала) от фактических значений забойного Р*з и значений забойного давления, вычисленных на основе модели Рз(а*(Д*)) на текущем интервале обработки КВД АЬ длительностью т, Ь0, - время
начала и завершения исследований; АҐ - интервал
обработки КВД, на котором показатель J принимает наименьшее значение.
Оценки параметров модели забойного давления «*(А* ) - потоков (1), (3), получены в (5) в момент времени f=Дf-тна основе метода адаптивной интерпретации КВД путем решения оптимизационных задач [6]
а,ф,) = ш^тпФ# ,$ ,), (6)
р: = argmin JG(a *p tD,
Pt
(7)
[Z> = Faa, + Г ,t tn -T]
(В)
на текущем интервале обработки Д^ Р*„ - век-
торы фактических забойных давлений и забойных давлений, вычисленных на основе модели КВД (1), (3), заданной с точностью до параметров «=(аи,аа); 1, векторы экспертных оценок
параметров пласта 1 и полученные на основе модели £>„; фНФ»Фа,Ю- оценка вектора управляющих параметров опр-еделяющих значимость (вес) экспертных оценок 1=(ги,га,-т); Fo, ¥а - матрицы. Например, для модели забойного давления радиального потока (1)
ДРз(,) = «1 +а21п(,), «1 = С[11'[К 1п\ 2,25к*? 1,
кь У т !лг )
CsquB
kxyL
(9)
и модели экспертных оценок вида
zit = ait = a1t + r1t, z 2t = an =a2t + r2,
— C quB f 2,25k ^ — C quB ait = ^ ln|------------4- |, aj = (10)
kxyL
kxyL
матрица ¥0 размерности (2хпк) и матрица ¥а размерности (2x2) имеют вид
„ = ( 11 .... 1 ^ „ = ( 10^
0 У Х1 Х2 Хпк) 2,пк ’ “ ^0 ^ 2,2’
где пк - число измерений забойного давления в текущем интервале обработки Дt=t+т (5); х^Ы^), 1=1, пк; кху, г2р - экспертные оценки радиальной проницаемости и приведенного радиуса скважины.
Момент начала радиального потока согласно (5) определяем по формуле
= Д/* -т’ (11)
где trt - время начала текущего интервала обработки, на котором показатель качества J (3) принимает наименьшее значение.
Следует отметить, что возможна и другая форма диагностического критерия (5) для определения времени начала радиального потока вида
tr = argmin
da*t(pt)
dt
t є[t0, tn -Т]-
(12)
где Ф(а„,Д) - комбинированный показатель качества интегрированной системы моделей КВД вида
[7]
[ К = р0а, +$, ’
За момент времени начало потока принимается t=tД при котором абсолютная величина производной оценки параметра модели КВД прини-
мает минимальное значение. Обоснованность использования критерия (12) подтверждается фактом стабилизации оценок КВД, полученных методом адаптивной интерпретации при выходе на соответствующий режим течения [6] (рис. 4, линия 3). Оценка радиальной проницаемости к*щ^ согласно (9) и (12) определяется по формуле
СдиВ
К ,t = -
t (p: ) l
при t = tr
(13)
Аналогично (13) в силу линейности по параметрам уравнений забойного давления (3) имеет место оценка горизонтальной проницаемости вида
,, (С аВ)2 и , ... ..
ку,, = 777)2^-г-ф-, при / = 1, , (14)
(77) лта2,(ф )
где t^ определяется по формуле (12).
Часто на КВД горизонтальных скважин наблюдается ранний и поздний радиальный потоки, представленные одним уравнением забойного давления (1) и слабо проявляющиеся на диагностическом графике (рис. 2). В этой связи целесообразным является использование диагностического критерия
t = arg max
d2alt (p*)
dt2
.t є^,tn -т]. (1Б)
Здесь за момент времени начало потока Г принимает то значение времени исследований, при котором абсолютная величина второй производной оценки параметра КВД «МР *) принимает максимальное значение, что фактически приводит к резкому «всплеску» абсолютной величины второй производной оценки а2(Р *) (рис. 6, 8, линия 3).
Результаты интерпретации КВД с диагностикой потоков
Результаты интерпретации КВД и диагностики потоков горизонтальных скважин № 1434 и 2031 нефтяного пласта месторождения Тюменской области приведены на рис. 3-6 и в табл. 1, 2. В качестве моделей КВД использовалось уравнение (1) и (3). Комбинированный показатель качества Фt интегрированной системы моделей КВД
= F,a, +$,,
|«2.' = a 21 +Г ,t tk -Т]
(1б)
был выбран в виде суммы частных квадратичных показателей качества
2
Ф, — (^ — ||АР3; - ^а, \\щ ) + р (а2,і - а2( )2, (17)
где экспертная оценка а2 і для радиального притока определялась из (10), а для линейного потока выбиралась равной
а2, —
4пЬк\ тк
Оценки радиальной к*у( и горизонтальной ку,( проницаемости рассчитывались по формулам
слуб
к' =■
к У, =
,(р) Ь' (СдБ)2 у
(Ьк) пта21(р{ )
при ,— ґг
при , =
(18)
где приближение параметра а% (А) моделей КВД (1) и (3) определялось из решения системы линейных алгебраических уравнений методом Гаусса [8]
(+ №) •(А ) = (FJWt + А а 2,0, (19)
в которой матрица формировалась согласно (10) для радиального потока ипри х1=л-1 для линейного потока; Ж(=^а§(^), £=1,ге*) - диагональная матрица весовых функций w(t) - определяет вес забойных давлений ДРз*^ в текущий момент времени t на участке обработки Дt=t+т с числом точек забойного давления ДР*г, Ь=1,ик (пк=5).
Оценки управляющего параметра определялись путем решения оптимизационной задачи (7) с использованием квадратичного показателя качества
3 = ||ДР3; - (А )||2)
методом золотого сечения [9].
На рис. 3 приведены фактические (линия 1) и восстановленные значения забойного давления (линии 2-4) на основе адаптивного метода наилучшего совмещения (АМНС), который следует из (19) при А=0, адаптивного метода наилучшего совмещения с регуляризацией (АМНСР), следующего из (18) при а2^=0 , и адаптивного метода интегрированных моделей (АМИМ) (19). Экспертные оценки радиальной и горизонтальной проницаемости для скважины № 1434 в начальный момент времени t=t0 выбиралисьравными кху=ку=100 мД, а для скважины № 2031 кх!1=ку=1000 мД с последующим их уточнением
СаиБ
—2------ при t > t0,
к у,, — ■
а 2, (р, )Ь (СдБ) 2 у (Ьк )2пта21 (р1 )
Из рис. 1 видно, что модель КВД, полученная АМИМ, более точно воспроизводит фактические значения забойного давления по сравнению с моделями, полученными АМНС и АМНСР.
На рис. 4-6 приведены оценки радиальной проницаемости (18) ее первой и второй производ-
ных, полученные АМНС, АМНСР и АМИМ соответственно.
150
145
I-
ГО
© 140
і
0)
5 135
га
Ч
а 130 і
ю
га
СО
1 I—і і і і і і і—і—і—і—і і і і і і і—і і і і і і—г-г-1—і і і
0,00 0,02 0,07 0,41 2,60 26,01 130,51
Время, ч
Рис. 3. Фактические (линия 1) и восстановленные значения забойного давления скважины № 1434
Рис. 4.
Оценки радиальной проницаемости скважины № 1434. 1) АМНС; 2) АМНС; 3) АМИМ
Всемя,ч
Рис. 5. Производные оценок радиальной проницаемости скважины № 1434 по критерию (12). 1) АМНС; 2) АМНС; 3) АМИМ
Из рис. 4 видно, что начиная с определенного момента времени наступает стабилизация оценок радиальной проницаемости и возникает вопрос об определении времени начала радиального притока для выбора подходящей оценки. Ответ на этот вопрос получен с использованием критериев (12) и (15), где первые и вторые производные радиальной проницаемости, изображенные на рис. 5 и 6 соответственно, дважды принимают минимальные и максимальные значения и фактически определяют момент времени начала раннего и позднего ра-
диальных потоков. Конкретные оптимальные значения времени начала этих притоков приведены в табл. 2.
0,00 0,02 0,07 0,41 2,60 26,01 130,51
Время, ч
Рис. 6. Вторые производные оценок радиальной проницаемости скважины № 1434 по критерию (15). 1) АМНС; 2) АМНС; 3) АМИМ
Следует отметить, что выделить ранний и поздний радиальные потоки по критериям (12) и (15) для скважины № 1434 при использовании АМНС и АМНСР не представлялось возможным.
На рис. 7, 8 для скважины № 2031 приведены абсолютные значения оценок первой (12) и второй (15) производных параметра а2»(А*) модели КВД (1), полученные АМНС, АМНСР и АМИМ соответственно.
0,01 0,08 0,39 1,70 7,25 30,64 129,25
Время,ч
Рис. 7. Производные оценок радиальной проницаемости скважины № 2031 по критерию (12). 1) АМНС;2) АМНС; 3) АМИМ
Время,ч
Рис. 8. Вторые производные оценок радиальной проницаемости скважины № 2031 по критерию (15). 1) АМНС; 2) АМНС; 3) АМИМ
Из рис. 7, 8 видно, что АМИМ с использованием критериев (12) и (15) позволяет диагностировать время начала раннего и позднего радиального
потоков. Конкретные оптимальные значения времени начала этих потоков приведены в табл. 2.
В табл. 1, 2 приведены результаты сравнительного анализа оценок радиальной и горизонтально проницаемостей и времени начала потоков, полученные при интерпретации горизонтальных скважин № 2031 и 1434 с использованием зарубежного программного комплекса 8арЫг, аналитического метода по формулам (2), (4) и АмНс, АМНСР и АМИМ.
Из табл. 1 видно, что оценки проницаемости нефтяного пласта при интерпретации КВД скважин № 2031 с использованием АМИМ с диагностикой потоков по критериям (12) и (15) согласуются с оценками, полученными с использованием программного комплекса Saphir. Для скважины № 1434 оценки проницаемости, полученные АМИМ, больше согласуются с АМНСР и АМИМ.
Преимущество АМИМ заключается в том, что он дает возможность определять время начала потоков в случае, когда использование аналитического (по формулам (2), (4)) и графоаналитического методов (рис. 2) вызывает значительные трудности.
Таблица 1. Результаты интерпретации скважин N 2G31 и 1434
Номер скважи- ны Метод интерпретации Оценки проницаемости, мД
Ранней радиальной Поздней радиальной Линей- ной
2031 Saphir - 3130 3742
АМНС 3161 3767 3341
АМНСР 2984 3087 3098
АМИМ 2849 3065 2803
1434 Saphir - 214 110
АМНС 842 1414 -
АМНСР - 1594 -
АМИМ 995 937 237
Таблица 2. Результаты диагностики радиального и линейного потоков скважин № 2031 и 1434
Номер скважины Метод диагностики потоков Время начала потоков, ч
Раннего радиального Позднего радиального Линей- ного
2031 Аналитический - 1,5 3,17
АМНС 0,21 1,7 2,28
АМНСР 0,21 2,28 3,05
АМИМ 0,11 7,25 9,58
1434 Аналитический - 0,75 1,89
АМНС 2,60 130,51 -
АМНСР - 51,90 -
АМИМ 0,65 32,75 51,90
Выводы
1. Предложен адаптивный метод диагностики радиального и линейного потоков с одновременной оценкой фильтрационных параметров нефтяного пласта в процессе гидродинамических исследований горизонтальных скважин по кривой восстановления забойного давления.
2. На примере интерпретации кривой восстановления давления двух горизонтальных сква-
жин однородно-пористого нефтяного пласта показано, что метод адаптивной интерпретации с диагностикой потоков дает возможность определять время начала радиального и линейного потоков в условиях, когда выделения потоков аналитическим способом и методом анализа
производной забойного давления вызывает значительные трудности.
3. Учет и корректировка экспертных оценок проницаемости нефтяного пласта позволяет повысить надежность определения времени начала радиального и линейного потоков горизонтальных скважин.
СПИСОК ЛИТЕРАТУРЫ
1. Шагиев Р.Г. Исследование скважин по КВД. - М.: Наука, 1998. - 304 с.
2. Иктисанов В. А. Определение фильтрационных параметров пластов и реологических свойств дисперсных систем при разработке нефтяных месторождений. - М.: ОАО «ВНИИОЭНГ», 2001. - 212 с.
3. Joshi S.D. Horizontal well technology. - Oklahoma: PenWell publ. comp., 1991. - 381 p.
4. Butler J.M. Horizontal wells for the recovery of oil, gas and bitumen. - Westmount: Petroleum Society Monograph, 1997. - 224 p.
5. Bourdet D., Ayoub J.A., Pirard Y.M. Use of pressure derivative in well test interpretation // SPE. - 1984. - № 12777. -293-302 p.
6. Гаврилов К.С., Сергеев В.Л. Адаптивная интерпретация нестационарных гидродинамических исследований скважин в системе «пласт-скважина» методом интегрированных моделей // Известия Томского политехнического университета. -2012. - Т. 321. - № 5. - С. 72-75.
7. Сергеев В.Л. Интегрированные системы идентификации. -Томск: Изд-во Томского политехнического университета, 2011. - 198 с.
8. Воеводин В.В. Вычислительные методы линейной алгебры. -М.: Наука, 1977. - 304 с.
9. Пантелеев А.В., Летова Т.А. Методы оптимизации в примерах и задачах. - М.: Высшая школа, 2002. - 544 с.
Поступила 03.06.2013 г.
UDC 519.688:622.276.5.001.42
ADAPTIVE INTERPRETATION OF PRESSURE RECOVERY CURVE OF HORIZONTAL WELLS WITH DIAGNOSTIC FLOW
E.V. Romanova, V.L. Sergeev Tomsk Polytechnic University
Relevance of the work is connected with the need for flow diagnostic at adaptive interpretation of unsteady hydrodynamic studies of horizontal wells using analytical models of pressure recovery curve.
The main aim of the study: to develop a method of adaptive interpretation of the pressure recovery curve which allows diagnosing streams at hydrodynamic studies of horizontal wells with simultaneous evaluation of reservoir filtration characteristics.
The methods used in the study: the theoretical and practical developments in the field of well testing, system analysis, system identification in light of additional prior information, optimization of functions and linear algebra were used. Solution of the problems of flon diagnosing, pressure recovery curve identification and interpretation has catried out on the basis of production data of downhole pressure in the well with the filter expert estimates of formation parameters using a computer program Saphir.
The results. The authors have developed the adaptive method for interpreting the pressure recovery curve with the diagnosis of radial and linear flow with simultaneous assessment of filtration parameters of the oil reservoir at hydrodynamic researches of horizontal wells. By the example of processing the results of flow testing by the pressure recovery curve of two horizontal wells of uniform-porous oil reservoir it was shown that the adaptive method allows increasing the reliability of determining the time of radial and linear flow starts.
Key words:
Interpretation, diagnostics, hydrodynamic tests, analysis of oil wells, a-priori information, horizontal wells, oil pools.
REFERENCES
1. Shagiev R.G. Issledovanie skvazhin po KVD [Well surveying by pressure build-up]. Moscow, Nauka, 1998. 304 p.
2. Iktisanov V.A. Opredelenie filtratsionnykh parametrov plastov i reologicheskikh svoystv dispersnykh system pri razrabotke neftya-nykh mestorozhdeniy [Determination of reservoir filtration parameters and dispersed system rheological properties when developing oil fields]. Moscow, BNIIOENG Publ., 2001. 212 p.
3. Joshi S.D. Horizontal well technology. Oklahoma, PenWell Publ.comp., 1991. 381 p.
4. Butler J.M. Horizontal wells for the recovery of oil, gas and bitumen. Westmount, Petroleum Society Monograph, 1997. 224 p.
5. Bourdet D., Ayoub J.A., Pirard Y.M. Use of pressure derivative in well test interpretation. SPE, 1984, no. 12777, 293-302 p.
6. Gavrilov K.S., Sergeev V.L. Adaptivnaya interpretatsiya nestat-sionarnykh gidrodinamicheskikh issledovaniy skvazhyn v siste-me «plast-skvazhyna» metodom integrirovannykh modeley [Adaptive interpretation of transient well test in the «layer-well» method of integrated models]. Bulletin of Tomsk Polytechnic University, 2012, vol. 321, no. 5, pp. 72-75.
7. Sergeev V.L. Integrirovannye sistemy identifikatsii [Integrated Systems of Identification]. Tomsk, TPU Publ. House, 2011. 198 p.
8. Voevodin V.V. Vychislitelnye metody lineynoy algebry [Computational methods of linear algebra]. Moscow, Nauka, 1977. 304 p.
9. Panteleev A.V., Letova T.A. Metody optimizatsii v primerakh i zadachakh [Optimization methods in the examples and problems]. Moscow, Vyshaya shkola, 2002. 544 p.
