CC BY
54
Вычислительные технологии
Том 11, № 3, 2006
ЧИСЛЕННАЯ ОПТИМИЗАЦИЯ РАЗМЕЩЕНИЯ СКВАЖИН
К. Р. АйДА-ЗАДЕ Азербайджанская государственная нефтяная академия, Баку e-mail: kamil_aydazade@rambler.ru
А. Г. БАГИРОВ Институт кибернетики НАН Азербайджана, Баку
Optimization problems for locations and extraction regimes of the newly commissioned wells at an oil field are considered. Mathematical formulation along with a numerical method for solving these problems is presented. The proposed numerical algorithm relies on the methods of optimal control over a system with distributed parameters.
Введение
Исследуются задачи оптимального размещения и управления дебитами скважин нефтяного пласта. Предполагается, что пласт является горизонтальным, однофазным, нефть — слабосжимаемой, тогда процесс фильтрации нефти можно принять двумерным [1, 2].
Эти задачи возникают при разработке новых и допроектировании эксплуатируемых месторождений, а также при оперативном управлении, планировании работы нефтепромысла. Математическая постановка описывается как параметрическая задача оптимального управления распределенными системами с сосредоточенными источниками, в которой параметрами служат координаты размещаемых скважин, управляющими воздействиями являются дебиты всех скважин нефтепромысла, ограничения и критерии оптимизации формируются с учетом технологических, геологических, экономических показателей эксплуатации месторождения. Применительно к разработке нефтегазовых месторождений подобные задачи исследовались многими авторами в различных вариантах технологических постановок, предложены и исследованы разные подходы численного решения [3-7].
В данной работе проведены исследования численного решения задачи оптимального размещения сосредоточенных источников — скважин (кустов скважин), получены выражения градиента функционала дискретизированной задачи по оптимизируемым параметрам, что позволяет использовать в дальнейшем для решения задачи эффективные методы оптимизации первого порядка. Приводится анализ численных результатов, полученных при решении контрольных задач.
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
1. Постановка задачи
Сформулируем задачу оптимального размещения скважин и оптимизации режимов их работы. Для описания процесса двумерной фильтрации слабосжимаемой нефти в пористой среде используем следующую краевую задачу [1, 2]:
с(х)|р - а1у(а(х)Ур) + ^ д1 (¿)5(х - х1) = 0, х € П° С Е2, * € (0; Т]; (1)
1=1
Р(х,^)|жег1 = Р1(х,*)
р(х, 0) = р°(х), х € П; (2)
ф(х, ¿)
= Р2(х,*), * € (0; Т], (3)
хеГ2
а(х) = к{х)к{х), ^(х) = + вп).
¡л
Здесь р = р(х, ¿) — давление в точке х € П в момент времени ¿; П — область фильтрации с границей Г, состоящей из непересекающихся частей Г1; Г2: Г = Г1 и Г2, Г1 П Г2 = 0 П° = П\Г; &(х) — проницаемость; ^(х) — мощность пласта; ^ — вязкость жидкости т — пористость; в], вп — коэффициенты сжимаемости жидкости и пористой среды х1 = (хг1;х2) — координаты мест размещения /-й скважины с дебитом д1 (¿), / = 1, 2, Ь — количество скважин; $(.) — обобщенная двумерная функция Дирака; Т — период времени планирования.
В рассматриваемой ниже задаче предполагается, что все функции, параметры пласта, нефти, участвующие в начально-краевой задаче (1)-(3), в том числе координаты скважин х1, / = 1,..., , заданы, а Ь2 скважин с неизвестными координатами хг, г = + 1,...,Ь1 + Ь2 = Ь, необходимо ввести в эксплуатацию, соблюдая следующие режимно-технологичес-кие, геологические, плановые показатели:
(х1,х2) € П, г = ¿1 + 1,...,Ь; (4)
||х" - х]|| > Д г,^ = 1,...,Ь, г = ; (5)
0 < д1 < д1 (*) < д1, / = 1,...,Ь; (6)
ь т
/ д1 > др1ап, (7)
1=10
где ||'|| — евклидова норма на плоскости; Д — минимально допустимое расстояние между скважинами; др1ап — плановое задание на добычу нефти. Как видно из (5)-(7), при решении задачи размещения новых скважин необходимо учитывать и оптимизировать также и режимы эксплуатируемых скважин, хотя может случиться, что дебиты некоторых работающих скважин заданы и менять их не рекомендуется.
В качестве критерия оптимальности в зависимости от производственной ситуации на промысле могут выступать такие показатели, как минимизация потерь пластовой энергии, максимизация добычи сырья и другие, включая их комбинации и многокритериальные случаи [5, 6, 8]. Рассмотрим задачу минимизации отклонения от средневзвешенной остаточной пластовой энергии:
/п £ т £
/ [р(х, Т) - р]2 ¿х + £1 ^ / [дг(*)] 2 ^+£2 X] ||хг||2, (8)
п 1=1 ° 1=£1+1
p = J J p(x,T )dx/mesH, n
где второе и третье слагаемые в (8) введены с целью регуляризации функционала; £1,£2 > 0 — параметры регуляризации.
Дополнительно к указанным выше условиям задачи приведем следующие важные замечания, исходящие из практических соображений и позволяющие существенно уменьшить объем вычислений при численном решении задачи. Предположим, что дебиты скважины — кусочно-постоянные во времени функции:
ql(t) = Qlk = const, t G [Tk-i,Tk), к =1,...,n, то = 0, r,n = T, l = 1,...,L, (9)
где n — число интервалов постоянства значений дебитов по всем скважинам.
Возможен вариант, когда каждая скважина (группа скважин, в частности куст скважин) имеет свои интервалы постоянства режимов:
ql(t) = qk, t G [rl-i,rl), rl0 = 0, т1щ = T, к =1,...,Ml, l = 1,...,L, (10)
где rо, ■■■,rlMl — времена переключения режимов для l-й скважины; Mi — число переключений. Нижеприводимые формулы в этом случае изменятся незначительно. Интервалы, определяемые значениями времен переключений режимов в (9), (10), могут быть заданы либо должны оптимизироваться. В данной работе будем требовать выполнения условия (9), а времена r1,...,rn — считать заданными. Практичность условия (9) следует из сложности реализации переменной во времени величины добычи сырья по скважинам, кустам скважин, что потребовало бы использования специальных телемеханических средств управления и регулирования, установленных на каждой скважине.
Размерность задачи оптимизации также существенно уменьшится, если агрегировать оптимизируемые параметры, а именно предположить, что ql(t) является дебитом l-го куста в целом, включающего несколько скважин (пять — десять).
Поставленная задача (1)-(9) является параметрической задачей оптимального управления системами с распределенными параметрами [3], в которой оптимизируются векторы:
Q = (Q1 ,...,Qn,...,QL,...,QL) G ELn, X = (xf1+1,xL1+1,...,xf,xL) G E2L2
общей размерности Ln + 2L2. Для ее численного решения предлагается использовать следующую схему. Применяя какую-либо конечно-разностную аппроксимацию всей задачи, необходимо получить конечномерную задачу математического программирования с ограничениями типа равенств и неравенств специального вида, которая относится к классу оптимизационных задач сетевой структуры [9, 10]. Для решения полученной задачи математического программирования следует использовать методы первого порядка, точнее, комбинацию метода внешних штрафных функций для учета ограничений (4), (5) и проекции сопряженного градиента штрафного функционала для учета линейных и позиционных ограничений (6), (7).
2. Конечномерная аппроксимация задачи
Введем равномерную сеточную аппроксимацию области П х [0, T], предполагая для простоты записи формул, что область П есть прямоугольник: [0, а] х [0, b], а Г2 = 0:
xu = ih1, x2j = jh,2, ts = sho, h = a/N1, h,2 = b/N2, ht = T/Nt,
и = {(i,j,s) : i = 0,...,Ni, j = 0,...,N2, s = 0,...,Nt} , причем тх = tSx, x = 1,..., n. Обозначим:
p(xii,X2j,ts)= Pijs, 9Z(ts) = QX = const, s = sx_b...,sx - 1, x = 1,...,n.
Пользуясь какой-либо схемой метода сеток [2, 11], аппроксимированную краевую задачу (1)-(3) запишем в достаточно общем виде:
Pjs = qijs), (i, j, s) G u; (11)
L
j = j (12) i=i
Зависимость (11), определенная функцией и ее множеством аргументов
^ijs = {Rnx : (f, П x) G uijs с u} , (13)
зависит от конкретной применяемой схемы метода сеток, а соответствующее индексное множество определяется выбранным шаблоном аппроксимации.
Формула (12) определяется схемой аппроксимации точечных источников. При этом, учитывая, что в общем случае их места расположения не совпадают с узлами сеточной области, их дебиты QS распределяются по близлежащим узлам с величинами, зависящими от расстояния до источника и обозначаемыми И, напротив, величины дебитов в узлах определяются суммой долей дебитов близлежащих источников, т. е. формулой (12). Ниже приведены некоторые возможные варианты схем аппроксимации точечных источников и их анализ.
Функция если для аппроксимации (1)-(3) используется неявная по t схема метода сеток, после несложных преобразований определятся следующим образом:
Pijs-1 + lPj+1s +
p..
L
+aijPij-is)/h2 - E qijs]/Aij, (i, j) G u0, s = 1,
1=1
(14)
, (i,j,s) G uij0,
I У (г,.) е г, 8 =
Ау = 1 + Л-о^г+у Л| + (Л1 + Л|)ау + ау+^Дсу Л^), а множество ее аргументов при (г,., в) € есть
^У« {_Рув—Ъ Ру^ рг— у«? Рг+У^ Ру — 1s, Ру+1«} • Таким образом, индексное множество ^у« в случае (14) есть
^ = {(г,.> - 1), ^.М^ (г - ^.М^ (г + (г. - М^ (г,. + М)} •
Для аппроксимации значений точечных источников (12) возможно использовать линейную функцию, функцию Гаусса, тригонометрическую (тангенсальную) функцию и др. [2, 3].
Пусть (иг,щ) определяют ячейку сеточной области, содержащую 1-ю скважину, т.е. имеют место неравенства
VIкг < х\ < (и + 1)къ ЩьЪ,2 < х12 < (щ + 1)к2.
Обозначим через V1 = {(иг,щг), (иг + 1,щ), (иг+ 1), (и+1,щ + 1)} множество номеров узловых точек сеточной области, ближайших к расположению 1-й скважины. Тогда линейную аппроксимацию дебита источника для Ь £ [тх-1,тх), в = вх-1,вх — 1, учитывая условие (9), можно представить в виде [2]
Я1х(х\ — Угк\)(х2 — щгк2)/(кг к2)2, г = иг + 1, Я1х((иг + 1)кг — х\)(х2 — щЬъ) / (к^)2, г Я1(х\ — игкг)((щ + 1)к2 — х\) / (к^)2, г
3 = Щг + иг, 3 = Щг + 1, иг + 1, 3 = Щг,
Я1х((иг + 1)кг — х\)((щ + 1)к2 — х\) / (к^)2, г = иг, 3 = Щ, 0, (г,3) /V1.
(15)
Несложно проверить, что
= XI ^ = Qlх, £ [тх-1,тх)-Аппроксимация источников с помощью функции Гаусса имеет вид [3]
я1-
QX ехр { — [(хи — х\)2 + (х2] — 4)2]/(2а2)} /(2па2), (г,3) £ Ш 0, (г,3) £Шг,
(16)
для в = вх-1, ■■■, вх — 1, где постоянная а2 характеризуется той эффективной окрестностью Ш1 1-й скважины — ячейками сеточной области вокруг расположения этой скважины, на которую распространяется влияние распределения дебита скважины. То есть при малых значениях а2 распределение мощности источника проводится для небольшого числа ячеек сеточной области, при больших значениях а2 аппроксимирующая функция имеет большую рассредоточенность по узловым точкам ш.
Ограничения (6), (7) и критерий (8) можно аппроксимировать следующим образом:
Яг < ^х < 11 , X =1,...,П, I =1,..;Ц
(17)
ь п
Т,^х(тх — тх-1) >
г=1 х=1
(18)
N1 N2 ь п ь
1(Х& р) = кгк2XX вгз [рРгN — р]2 + е^Х )2(т* — т—)+ ¿2 ^ ||хг ||2, (19)
г=о j=o г=1 8=1 г=Ь!+1
N1 N2
р^/(N1X2).
г=о j=o
Здесь — коэффициенты применяемой квадратурной формулы.
Таким образом, получена задача конечномерного математического программирования (4), (5) (11), (17)-(19), в которой основными являются ограничения типа равенств (11).
Пользуясь структурой этих ограничений, запишем формулы градиента функции (19) по Ьп + 2/2 оптимизируемым параметрам:
V/(Х,д; р) = (Уд/(Х,д; р), Ух / (Х,д; р)). (20)
Введем вспомогательный вектор импульсов [9, 10]:
Щв = -^-, € ^ (21)
Здесь производные понимаются как полные с учетом того, что изменение значения давления руз в точке (ж1г, ) в момент времени ^ согласно зависимостям (11) повлечет за собой изменение давления в других точках Введем понятие сопряженного к индексного множества:
= : (г,з,в) € , (г,3,в) € и. (22)
Из (11), (12) и этого определения следует, что если (£,п,х) € то в (11) имеется зависимость р^ = (...,р^пх,...), а из € следует, что среди соотношений (11)
имеется зависимость р^пх = /£пх(...,руз,...). Тогда, учитывая (11), (12), соотношения (21) можно записать в следующем виде:
= д/(Х,Ф р) + у^ ^ р) др^х = = д/(Х,Ф р) + ^ д^пх.
+ Е ^р«пх, (г,з,в) € ^ (23)
= <
Например, для неявной схемы метода сеток (14), пользуясь конкретизированными для этого случая сопряженными индексными множествами (22), в частности для последнего временного слоя в = будем иметь
{ (г - 1,з, ЛТ4), (г + 1,з, ЛТ4), (г,;/ - 1, ЛТ4), (г,з + 1, ЛТ4), 2 < г < N1 - 2, 2 < з < N - 2, (1, 2, N4), (2,1,^), г = 1, з = 1,
(Ж1 - 2,1, N4), (^ - 1, 2, N4), г = N - 1, з = 1, (24)
(1,^ - 2, N4), (2, N2 - 1, N4), г =1, з = N2 - 1,
(^ - 2, N - 1, N4), (^ - 1, N - 2, N4), г = N - 1, з = N - 1,
0, (г,з) € Г,
а частная производная функционала (19), участвующая в (23), определяется непосредственным дифференцированием:
{0, в < N4, 0 < г < 0 < з < N¡2,
/ N1 N2 \
2^1 Ц ву Ьу^ - р] - Е Е вгкЬгВД - р]/^2) , (25)
\ г=0к=0 /
0 < г < N1, 0 < з < N2, в = N4.
Все частные производные вида др^пх/дргу5 в (22) непосредственно и несложно определяются из (14).
Пользуясь (14), (23)-(25), приведем соотношения для вспомогательных переменных в случае использования неявной схемы аппроксимации (13), например для s = Nt, i = 2,Ni - 2, j = 2,...,N2 - 2:
. _dl(X, Q; p) t дрг -ijNt , , dpi+ijNt . .
j = dpijNt Фi-ijNt + Фг+1т +
. dpij-iNt , . dpij+iNt ,
дрг-Ху81дргу8 = /к\, дрг+1^8/др^8 = Л/к\,
дрч-и/дргэз = Лч-1/Н22, д'[Щ+1з/д'[Щ3 =
Л%3 = а^/[1 + ко(а^ Ьк2 + а^ (к\ + Ьк2) + а^+хЬ^ )/(с^ к\к22).
Таким образом, в случае неявной аппроксимации (14) исходной задачи после определения значений из линейных соотношений, определенных формулой (23), находятся значения импульсов, но в обратном порядке по в: сначала для в = N1 решается система N1N2 уравнений относительно р^Яи я = 1, N1, ] = 1,..., N2, а потом — для в = N1 — 1 и т. д. до в = 0.
Учитывая, что дебиты принимают кусочно-постоянные значения, т.е. условие (9), и то, что работа скважины непосредственно (в явном виде) влияет на значения давления из (14) лишь в тех узлах сетки, которые участвуют в (15) или (16) (в зависимости от схемы аппроксимации функций источников), компоненты градиента функционала по управляющим параметрам — дебитам скважин и координатам их размещения — определятся следующими формулами:
-I д1 V- ^т^ дрц3 , , ,
-ох = ж +£ X , х = 1,...,п, 1 = 1,...,ц (2б)
-Г = Г Г ' 1 = 1х + к =1,2. (27)
к 3=0 1
Частные производные, участвующие в правых частях (26), (27), несложно определяются непосредственным дифференцированием из формулы (19) — сд1 /д01х, из (14), (15) —
др^з/д01х, др^.3/дх1.
Таким образом, значения компонентов градиента определяются по формулам (26), (27) после решения краевой задачи по формулам (14) и нахождения значений вспомогательных импульсных переменных по формулам (23). В целом для решения редуцированной конечномерной задачи оптимизации (11), (17)-(20) можно использовать методы первого порядка, о которых говорилось выше. Отметим, что приводимые выше формулы для градиента функционала могут быть получены методом Лагранжа. В этом случае импульсные переменные выступают в роли множителей Лагранжа.
3. Результаты численных экспериментов
Ниже приводятся результаты проведенных экспериментов по численному решению задач оптимизации размещения скважин и управления их дебитами. Численный анализ прово-
дился не только непосредственно над оптимизационными задачами, а исследовались вспомогательные вычислительные задачи, с которыми приходится сталкиваться при решении основной проблемы, и они существенно влияют на результаты решения всей задачи.
Приводимые численные результаты получены в основном при следующих условиях и значениях параметров задачи:
П==[0;3 • 105 (см)] х [0; 6 • 105 (см)], dp(x,t)/dn = 0, x 6 Г,
h(x) = const = 500 см, p(x, 0) = p0 (x) = const = 150 атм, в = 0.000225, вп = 0.000010, mo = 0.2, ^ = 2.5 спз,
k = 0.4 Дарси.
Все расчеты проводились с двойной точностью на алгоритмическом языке Паскаль на компьютере Pentium-1 с процессором 160 МГц.
Прежде всего сделаем важное для поставленной задачи замечание. Попытка использовать для задачи размещения скважин линейную аппроксимацию функции источников (14) не увенчалась успехом по той причине, что целевая функция I(X,Q;p) по первому аргументу во всех точках прямых, образующих сеточную область, недифференцируема.
На рис. 1, а приведен график целевой функции для тестовой задачи с одной скважиной, координаты которой изменялись по прямой x2 = 3 (км), а на рис. 1, б приведены линии уровня целевой функции в первой четверти декартовой системы, центр которой совпадает с центром прямоугольника — П, а оси параллельны сторонам. Из обоих рисунков видно сложное (невыпуклое, негладкое) поведение функционала по первому аргументу — координатам скважин, в то же время по второму аргументу — дебитам скважин — целевая функция ведет себя гладко. Такое поведение функции, как это видно из (14), (15), объясняется скачкообразным поведением аппроксимации при переходе от одной клетки сеточной области к другой.
Эти же функции приводятся на рис. 2 при использовании аппроксимации (16), когда координаты скважины изменяются, например, по прямой xi + x2 = 3 в первом ортанте,
Рис. 1.
Рис. 2.
начиная от центра. Эффективная окрестность Ш1 в этих расчетах состояла из девяти клеток и включала кроме V1 восемь клеток из окрестности V1. Величина а2 рассчитывалась в зависимости от к1,к2, а именно выбиралась из условия, что
1,5^1 1,5^2
1
2па2
exp[-(X + x^2)/(2a2)]dx1dx2 > 1 — £0,
- 1,5hi - 1,5h2
где величина > 0 связана с точностью проводимых расчетов.
В таблице приводятся результаты сравнения по точности и времени расчета (на Pentium 1) производных целевого функционала по координатам расположения пяти скважин, полученных с использованием формул (26), (27) и с помощью конечно-разностных аппрок-
симации:
(28) (29)
61(Я;р)/62г = I(Я + Ье{;р) - I(Я - 5вг;р)/(25) + 0(52);
61 (Я;р)/6Яг = I(Я + 5ег;р) - I(Я;р)/5 + 0(5).
Здесь ег — ¿-орта, г =1, 2,..., 2Ь, Ь = 5.
Отметим, что для подбора "хорошего" значения 5, используемого в формулах (28), (29), пришлось провести специальные эксперименты. В качестве параметров регуляризации в критерии (8) использовались значения е\ = 0.5 • 10-3-к, = 0.5 • 10-5-к, к = 0,1, 2, т.е. трижды проводилась последовательная оптимизация. Отметим, что регуляризация позволила получать устойчивое оптимальное значение при различных начальных значениях оптимизируемых параметров.
1 Координаты (dI/dx[,dI/dxl2) ■ 105
скважины (x1,4) Формула (27) Формула (28) Формула (29)
1 (2.2; 2.04) (0.59; —0.15) (0.64; —0.18) (0.66; —0.24)
2 (1.02;2.04) (0.158; —0.346) (0.16; —0.367) (0.164; —0.29)
3 (0.78;1.56) (—1.15; 0.418) (—1.172; 0.307) (—1.02; 0.27)
4 (1.23; 2.46) (—0.19; —1.10) (—0.17; —1.03) (—0.15; —0.99)
5 (1.83; 3.66) (0.08; —1.9) (0.073; —1.85) (0.11; —2.37)
Время счета, с 15 183 82
Рис. 3.
На рис. 3 приведены результаты оптимизации размещения скважин для случая шести скважин и различных начальных размещений (на рисунках начальное размещение скважины отмечено знаком " Ж", полученное оптимальное размещение — знаком " ■"). Относительно значения дебитов проводились разные эксперименты для различного числа интервалов постоянства их значений. В частности, для случая пяти скважин и п = 1 (т. е. все скважины во времени не меняют режим работы) получено оптимальное значение = 33, 250, г = 1,6.
Заключение
Исследуемый в работе подход к решению задач размещения скважин нефтепромысла и оптимального управления ими может быть использован в решении многих других проблем управления распределенными системами, в которых управляющими воздействиями являются сосредоточенные источники, в том числе и для случаев, если они перемещаются. Ясно, что используемая математическая модель, а именно в виде (1)-(3), в полученных формулах в принципе не имела большого значения, точно так же, как и применяемые конкретные схемы конечномерной аппроксимации. Преимущество данного подхода, полученных формул в том, что они позволяют использовать для решения задач богатый арсенал численных методов конечномерной оптимизации.
Список литературы
[1] ЧАРНЫй И.А. Подземная гидрогазодинамика. М.: ГНТИ, 1963.
[2] Азиз Х., Сеттари Э. Математическое моделирование пластовых систем. М.: Недра, 1982.
[3] Бутковский А.Г. Методы управления системами с распределенными параметрами. М.: Наука, 1975.
[4] Управление разработкой нефтяных месторождений / Я.М. Берщанский, В.Н. Кулибанов, М.В. Мееров, О.Ю. Пермин. М.: Недра, 1983.
[5] АхмЕтзянов А.В., Куливлнов В.Н. К проблеме оптимального управления нефтяных месторождений // Автоматика и телемеханика. 1999. № 6. С. 5-13.
[6] Ахметзянов А.В., КуливАнов В.Н. Оптимальное размещение источников для стационарных скалярных полей // Автоматика и телемеханика. 1999. № 6. С. 50-58.
[7] äyda-zade k.r., Handzel a.v. Optimization of oil field wells' placement // III Intern. Conf. "Math. and Comp. Appl." (ICMCA-2002). Turkey, 2002. P. 108-110.
[8] Авдуллаев Ф.М., Айда-заде К.Р., Калаушин М.А. Система управления Ямбургским промыслом // Приборы и системы управления. 1990. № 2. С. 32-37.
[9] АйДА-ЗАДЕ К.Р. Исследование и численное решение конечно-разностных аппроксимаций задач управления распределенными системами // Журн. вычислительной математики и мат. физики. 1989. № 3. С. 122-130.
[10] АйДА-ЗАДЕ К.Р. Исследование нелинейных оптимизационных задач сетевой структуры // Автоматика и телемеханика. 1989. № 4. С. 63-71.
[11] Самарский А.А. Теория разностных схем. М.: Наука, 1977.
Поступила в редакцию 3 мая 2005 г.
