Численное решение задач идентификации коэффициента фильтрации на основе двухшаговых методов минимизации функции невязки

При решении многих важных практических задач управления, идентификации и адаптации возникает нелинейная задача о наименьших квадратах где J (x) — функция невязки (целевая функция), х = (x,., xn) r — вектор переменных минимизации, R = (rx,., rm) T — вектор невязки, rt = ri (х) — нелинейные функции, п — число переменных минимизации, т — число невязок. При численном решении нелинейной задачи о наименьших квадратах, как правило, строразуют убывающую последовательность, т. е. значения целевой функции в точках последовательности удовлетворяют неравенству (условие убывания функции невязки) минимизации различаются способами построения релаксационной последовательности и в зависимости от порядка используемых производных делятся на три группы [3,4,7,44,47,51,56]: методы нулевого порядка (или методы прямого поиска), методы первого (градиентные методы) и второго порядка. Методы нулевого порядка используют только значения целевой функции. К методам нулевого порядка относятся метод конфигураций, метод деформируемого многогранника, метод Розенброка, методы случайного поиска и др. В методах первого порядка используются значения целевой функции и значения её частных производных первого порядка. Наиболее часто используемые методы этого типа: градиентный метод, метод наискорейшего спуска и метод сопряженных градиентов. В методах второго порядка (метод Ньютона, метод Ньютона с по.

14,19,90]: ится последовательность точек хА в которых значения целевой функции об.

J{xk)< J (xk~l).

Такая последовательность точек |хЧ называется релаксационной [3]. Методы иском шага) используются значения целевой функции и значения её частных производных до второго порядка включительно. При решении практических задач построение матрицы вторых производных (матрица Гессе) требует больших вычислительных затрат. Поэтому на практике обычно применяются методы квазиньютоновского типа. В этих методах используется аппроксимация матрицы вторых производных или её обратной с помощью частных производных первого порядка. При приближении к точке минимума эти методы по своим свойствам приближаются к методам второго порядка. К методам, в которых аппроксимируется матрица обратная к матрице Гессе, относятся методы Дави-дона-Флетчера-Пауэлла, Бройдена-Флетчера-Голдфарба-Шэнно, Пауэлла и др., в методе Гаусса-Ньютона аппроксимируется матрица Гессе.

Метод Гаусса-Ньютона лежит в основе многих методов квазиньютоновского типа и создан специально для минимизации функций, имеющих вид суммы квадратов. В этом методе матрица Гессе аппроксимируется’с помощью матдг.

— > - матрица чувствительности, и переход на новые cbcfJ значения переменных осуществляется, по формуле x" ewxo! d + dGN, dGN = -H~lg — направление спуска Гаусса-Ньютона, g — градиент функции невязки. В процессе-минимизации методом Гаусса-Ньютона матрица Н может оказаться плохо обусловленной или вырожденной, что вносит большие погрешности в определение направления спуска, и значение функции невязки при переходе на новые значения переменных может вырасти. Различные модификации метода Гаусса-Ньютона направлены на устранение этих недостатков (метод Гаусса-Ньютона с поиском шага, метод Левенберга-Марквардта).

В методе Гаусса-Ньютона с поиском шага используются различные процедуры поиска шага вдоль направления спуска. Длина шага определяется из условия минимума функции вдоль направления спуска с использованием методов одномерной минимизации или с помощью процедуры дробления шага. рицы Н = А А, где, А =.

Другой модификацией метода Гаусса-Ньютона является метод Левенбер-га-Марквардта. В этом методе последовательные приближения определяются формулой x" cw = xold — (Н + juE)~] g, где ju> 0 — параметр Марквардта. Эта формула была независимо предложена Левенбергом [78] и Марквардтом [81]. При больших значениях параметра Марквардта направление минимизации метода Левенберга-Марквардта близко к направлению метода наискорейшего спуска, при малых значениях — к направлению метода Гаусса-Ньютона. В начале процесса минимизации, когда текущая точка расположена вдали от точки минимума, выбираются большие значения параметра Марквардта. При приближении к точке минимума значения параметра Марквардта уменьшаются. Таким образом, метод Левенберга-Марквардта, в отличие от метода Гаусса-Ньютона, не требует хорошего начального приближения. Направление спуска метода Левенберга-Марквардта d ~-{Н + juE)~l g является решением задачи на условный минимум [14,19] mmi|| Ad + R2 при ограничении d2 < Д, где, А — матрица чувствительности, R — вектор невязок, А — параметр, связанный с параметром jli, определяет размер доверительной области. Поэтому метод Левенберга-Марквардта можно отнести к классу методов доверительной области [19]. Существуют различные модификации метода Левенберга-Марквардта, отличающиеся стратегиями выбора параметра Марквардта. В [44,90] начальное значение параметра Марквардта выбирается достаточно большим. В случае выполнения условия убывания функции невязки на итерации при переходе на следующую итерацию параметр Марквардта уменьшается, а при нарушении этого условия параметр Марквардта увеличивается до тех пор, пока оно не выполнится. В [44] коэффициент изменения параметра Марквардта берётся равным двум, а в [90] - десяти. В' [75] в начале каждой итерации параметр Марквардта равен нулю. В случае нарушения условия убывания функции невязки параметр Марквардта увеличивается по формуле /л = 1.5// + 0.001 до тех пор, пока это условие не выполнится. Этот вариант метода Левенберга-Марквардта используется в пакетах UCODE и MODFLOWP, предназначенных для решения обратных задач. Вариант метода Левенберга-Марквардта, в котором подбирается А, а // вычисляется по значению, А приводится в [86]. Этот вариант включен в пакет программ MINPACK. Процедура подбора значения параметра Марквардта в зависимости от отношения между фактическим изменением функции невязки в результате пробных шагов и ожидаемыми величинами этих изменений построена в [68]. В [34] параметр Марквардта определяется методом золотого сечения из условия минимума функции невязки.

Для проверки эффективности методов минимизации обычно используются специальные тестовые функции. Различные примеры тестовых функций можно найти в [7,59,60,87].

К задаче минимизации функции невязки сводится задача идентификации коэффициента фильтрации водоносных и нефтяных пластов при фильтрации в них жидкости по замерам напора в наблюдательных точках [75,90,92]. В этом случае функция невязки имеет вид суммы квадратов разностей между вычисленными и измеренными значениями напора в наблюдательных точках. Коэффициент фильтрации ищется в виде кусочно-постоянной функции. Область решения задачи разбивается на зоны однородности, каждая из которых характеризуется своими значениями коэффициента фильтрации [64,66,89]. Задача определения поля напора при известном коэффициенте фильтрации является прямой задачей, а задача идентификации коэффициента фильтрации по замерам напора в наблюдательных точках — обратной коэффициентной задачей.

Задача идентификации коэффициента фильтрации относится к классу некорректно поставленных задач [1,5,6,10,17,18,41,50,52,53,58,90]. Задача называется некорректно поставленной, если нарушается одно из условий: решение существует, решение единственно, решение устойчиво относительно изменения исходных данных. В задачах определения параметров водоносных пластов может нарушаться любое из условий корректности. Примеры таких задач приведены в [62,83,90,91]. При решении обратных коэффициентных задач из-за наличия погрешностей обычно на первых итерациях значения параметров приближаются, а, начиная с некоторых итераций, удаляются от своих истинных значений, при этом функция невязки продолжает уменьшаться. Для выбора номера итерации с более достоверными значениями коэффициента фильтрации применяются специальные правила останова процесса минимизации или прерывания полученной последовательности значений функции невязки. Эти правила являются одним из регуляризирующих элементов решения обратных задач [1,5,6,41,52].

При решении задач идентификации коэффициента фильтрации с большим числом идентифицируемых параметров обычно используются квазиньютоновские методы минимизации функции невязки. Для вычисления. градиента функции невязки и матрицы чувствительности используется один из методов: метод конечно-разностных соотношений, вариационный метод, метод прямого дифференцирования [90,92]. Недостатком метода конечно-разностных соотношений является, сложность вычисления производных с требуемой точностью. Более точные значения производных получаются при — использовании метода-прямого дифференцирования и вариационного метода. В [54] приведена общая схема построения вариационного метода, в [76] вариационный метод впервые использован при решении задач идентификации параметров. Сравнение трёх перечисленных выше методов вычисления производных при решении задачи идентификации коэффициента фильтрации приводится в [22]. Для нахождения^ матрицы чувствительности методом конечно разностных соотношений или методом прямого дифференцирования необходимо решить п+1 уравнений, аналогичных уравнению фильтрации (п — число идентифицируемых параметров), а в случае использования вариационного метода т+1 уравнений (т — число замеров напора).

В процессе решения задачи идентификации коэффициента фильтрации на каждой итерации требуется неоднократно решать прямую задачу. Решение прямой задачи проводится численными методами. Система дифференциальных уравнений с частными производными (уравнение фильтрации и соответствующие начальные и граничные условия) сводится к системе линейных алгебраических уравнений. Для этого используются метод конечных разностей, метод конечных элементов, метод контрольных объёмов и др. [8,15,26,36,40,42,4850,55,65,69,71,72,79]. В практических задачах получаемая система алгебраических уравнений имеет большую размерность и обычно решается одним из итерационных методов, учитывающих симметричность и сильную разреженность матрицы системы. Одним из наиболее часто используемых методов является метод сопряженных градиентов с предобусловливанием [74,77,84,85]. Скорость сходимости и эффективность этого метода сильно зависят от используемой предобусловливающей матрицы. Сравнение результатов, полученных методом сопряженных градиентов с использованием трёх различных предобусловли-вающих матриц (неполное разложение Холесского, модифицированное неполное разложение Холесского, полиномиальный метод), проведено в [74]. В [77] приведены результаты решения девяти трёхмерных задач по определению поля напора в условиях однофазной стационарной фильтрации жидкости, полученные методом сопряженных градиентов без предобусловливания и с использованием для построения предобусловливающей' матрицы диагонального масштабирования, неполного разложения Холесского, неполной факторизации, модифицированной неполной факторизации. Другие подходы к построению предобусловливающей матрицы используются в [70,80,88].

Для тестирования методов решения задачи идентификации коэффициента фильтрации используются модельные задачи [63]. В модельных задачах по известным значениям коэффициента фильтрации (истинные значения) определяются значения напора в наблюдательных точках. Далее значения коэффициента фильтрации считаются неизвестными, и требуется их определить по значениям напора в наблюдательных точках. Число и границы зон однородности могут считаться как известными, так и неизвестными. В последнем случае их также требуется определить в процессе идентификации. Методы решения отрабатываются на модельных задачах, как без погрешностей, так и с погрешностями, вносимыми в задачу. В модельных задачах без погрешностей решение всегда существует. Число и положение наблюдательных точек должно быть достаточным для выполнения условия единственности решения в задаче без погрешностей, что проверяется численными экспериментами. Полученные врезультате решения модельной задачи идентификации значения коэффициента фильтрации всегда можно сравнить с истинными значениями, и по результатам сравнения оценить эффективность метода решения задачи. При решении практической задачи идентификации коэффициента фильтрации приходится многократно строить и решать различные модельные задачи. При решении модельных задач приходится определять достаточность исходных данных для идентификации того или иного набора параметров, влияние погрешностей в исходных данных и используемых регуляризирующих процедур. Всё это требует больших вычислительных затрат.

Основной целью данной работы является разработка методов*минимизации функции невязки, сокращающих вычислительных затраты. В работе предложены двухшаговые методы минимизации функции невязки: двухшаговые методы Ньютона, двухшаговые методы Гаусса-Ньютона и двухшаговые методы Левенберга-Марквардта. В двухшаговых методах первый <, шаг каждой итерации проводится по алгоритмам классических методов, но допускается-увеличение функции невязки. Итоговые же значения функции невязки на итерациях, как и в классических методах, образуют убывающую последовательность. Для анализа известных методов минимизации функции невязки и построения двухшаговых методов введено понятие запаса чувствительности. При построении двухшаговых методов используется главная система координат в пространстве переменных минимизации. Главная система координат вводится с помощью сингулярного разложения симметрической приближённой матрицы вторых производных [16,19]. На первом этапе сингулярного разложения с помощью преобразований Хаусхолдера строится подобная трёхдиагональная симметрическая матрица, которая на втором этапе итерационно приводится к диагональному виду с помощью QR-итераций с одинарным сдвигом. В главной системе координат запас чувствительности характеризует потенциальную возможность к минимизации функции невязки вдоль каждого из направлений.

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

В [54] приводится пример решения системы нелинейных уравнений с двумя неизвестными методом Ньютона. В этом примере на каждой итерации проводится масштабирование невязок, что позволяет сократить числом итераций, необходимое для получения решения. На каждой итерации процесса* минимизации значение функция невязки с масштабными коэффициентами уменьшается, однако, значение функция невязки без масштабных коэффициентов на отдельных итерациях увеличивается.

Для решения более сложных задач (большое число переменных минимизации, плохо обусловленные задачи) обычно используется метод Левенберга-Марквардта. На основе различных вариантов метода Левенберга-Марквардта построены двухшаговые методы Левенберга-Марквардта. В двухшаговых методах Левенберга-Марквардта используется главная система координат, полученная с помощью сингулярного разложения приближенной матрицы вторых производных. Направления минимизации в главной системе координат условно делятся на две группы: направления, соответствующие большим сингулярным числам, и направления, соответствующие маленьким сингулярным числам. Смещение по направлениям первой группы соответствует спуску ко дну оврага, по направлениям второй группы — движению вдоль склона оврага. На первом шаге двухшаговых методов Левенберга-Марквардта допускается увеличение функции невязки за счёт её увеличения в направлениях с большими сингулярными числами. Вдоль этих направлений проводится смещение параметров на втором шаге. В модифицированных двухшаговых методах Левенберга-Марквардта проводятся дополнительные смещения переменных вдоль направлений с большими сингулярными числами.

Предложенные двухшаговые методы минимизации функции невязки относятся к монотонным методам. В работах [67,93] рассмотрены немонотонные методы, в которых допускается увеличение функции невязки на отдельных итерациях. В немонотонных методах при переходе на новую итерацию значение функции невязки сравнивается не с предыдущим значением, а с максимальным значением из М предыдущих итераций, где М — задаваемый параметр.

Численные алгоритмы, основанные на двухшаговых методах Левенберга-Марквардта, протестированы на шести модельных задачах идентификации коэффициента фильтрации трехмерного анизотропного пласта по замерам напора в наблюдательных точках в условиях однофазной стационарной фильтрации жидкости, подчиняющейся закону Дарси. Модельные задачи представляют собой водоносный пласт, состоящий из пяти слоев. Каждый слой разбит на зоны однородности коэффициента фильтрации (71 зона однородности). Каждая зона однородности характеризуется двумя значениями коэффициента фильтрации. В задачах требовалось восстановить 142 значения коэффициента фильтрации. На кровле пласта заданы граничные условия второго рода. Подошва и боковая поверхность пласта считались непроницаемыми, за исключением участков боковой поверхности пятого слоя. Модельные задачи отличаются числом и положением наблюдательных точек, граничными условиями. Для определения поля напора уравнение фильтрации с соответствующими граничными условиями методом конечных элементов Галёркина аппроксимировалось системой линейных алгебраических уравнений. Полученная система уравнений решалась методом сопряжённых градиентов с предобусловливанием в виде неполного разложения Холесского. Значения коэффициента фильтрации определялись из минимума функции невязки, имеющей вид суммы квадратов разности между вычисленными и измеренными значениями напора в наблюдательных точках. Сравнение результатов решения модельных задач идентификации коэффициента фильтрации классическими и двухшаговыми методами показало^ чтопо вычислительным затратам наиболее эффективными являются двухшаговые методы Левенберга-Марквардта с дополнительными спусками ко дну оврага. При решении задач с погрешностями в замерах напора значения" идентифицируемых параметров, как правило, начиная с некоторой итерации, удаляются от своих истинных значений, хотя функция невязки продолжает уменьшаться. Для получения итоговых значений коэффициента фильтрации более близких к истинным, в работе использовался критерий выбора номера итерации с итоговыми значениями идентифицируемых параметров.

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

Цели диссертационной работы:

— построение эффективных методов минимизации функции невязки;

— разработка на основе предложенных методов минимизации функции невязки численных алгоритмов решения задачи идентификации коэффициента фильтрации водоносного пласта.

Структура и краткое содержание работы.

Диссертация состоит из введения, четырёх разделов, заключения и списка литературы. Работа содержит: страниц — 150, рисунков — 46, таблиц — 26, список литературы — 93 наименование.

Выводы.

1. На основе методов Левенберга-Марквардта, Гаусса-Ньютона и двухшагового метода Левенберга-Марквардта построены методы минимизации функции невязки, учитывающие сравнительную априорную информацию о значениях идентифицируемых параметров. Особенностью этих методов является то, что упорядоченность параметров сохраняется в течение всего процесса минимизации. Число переменных минимизации может быть меньше числа идентифицируемых параметров и увеличиваться постепенно в процессе минимизации.

2. Методы минимизации функции невязки, учитывающие сравнительную априорную информацию о значениях идентифицируемых параметров, протестированы при численном решении модельных задач идентификации коэффициента фильтрации. Учёт априорной сравнительной информации о значениях идентифицируемых параметров позволил сократить вычислительные затраты и получить итоговые значения коэффициента фильтрации более близкие к истинным. Из методов, учитывающих априорную сравнительную информацию, наиболее эффективным по вычислительным затратам показал себя двухшаговый метод.

Заключение

.

В диссертации построены двухшаговые методы минимизации функции невязки: двухшаговые методы Ньютона, двухшаговые методы Гаусса-Ньютона и двухшаговые методы Левенберга-Марквардта. В двухшаговых методах первый шаг каждой итерации проводится по алгоритмам классических методов, но допускается увеличение функции невязки. Итоговые же значения функции невязки на итерациях, как и в классических методах, образуют убывающую последовательность.

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

В двухшаговых методах Левенберга-Марквардта используется главная система координат, полученная с помощью сингулярного разложения приближенной матрицы вторых производных. Направления минимизации в главной системе координат условно делятся на две группы: направления, соответствующие большим сингулярным числам, и направления, соответствующие маленьким сингулярным числам. Смещение по направлениям первой группы соответствует спуску ко дну оврага, по направлениям второй группы — движению вдоль склона оврага. На первом шаге двухшаговых методов Левенберга-Марквардта допускается увеличение функции невязки за счёт её увеличения в направлениях с большими сингулярными числами. Вдоль этих направлений проводится смещение параметров на втором шаге. В модифицированных двухшаговых методах Левенберга-Марквардта проводятся дополнительные смещения переменных вдоль направлений с большими сингулярными числами.

На примерах минимизации тестовых функций (квадратичная функция простой структуры, функция Розенброка, двумерная экспоненциальная функция, функции Биля, Пауэлла и Зангвилла) проведено сравнение известных методов: Ньютона, Гаусса-Ньютона, Ньютона и Гаусса-Ньютона с поиском шага, Давидона-Флетчера-Пауэлла, Бройдена-Флетчера-Голдфарба-Шэнно, Левенберга-Марквардта с двухшаговыми методами Ньютона и Гаусса-Ньютона без поиска и с поиском шага. Результаты тестирования показали, что двухшаговый метод Гаусса-Ньютона с поиском шага позволил достичь точки минимума для всех рассмотренных тестовых функций с наименьшими вычислительными затратами.

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

На основе методов Левенберга-Марквардта, Гаусса-Ньютона и двухшагового метода Левенберга-Марквардта построены алгоритмы минимизации функции невязки, учитывающие априорную сравнительную информацию о значениях идентифицируемых параметров. Учёт априорной сравнительной информации о значениях идентифицируемых параметров позволил сократить вычислительные затраты и получить итоговые значения коэффициента фильтрации более близкие к истинным. Из этих методов наиболее эффективным по вычислительным затратам показал себя двухшаговый метод.