CC BY
24
УДК 519.653:519.63
DOI: 10.18698/1812-3368-2016-6-42-55
ПРИМЕНЕНИЕ БЫСТРОГО АВТОМАТИЧЕСКОГО ДИФФЕРЕНЦИРОВАНИЯ ПРИ НАХОЖДЕНИИ ИСПАРЕНИЯ С ПОВЕРХНОСТИ ПОЧВЫ
B.В. Дикусар1 dikussar@yandex.ru
C.В. Засухин2 s.zasukhin@yandex.ru
1 Вычислительный центр им. А.А. Дородницына ФИЦ «Информатика и управление» РАН, Москва, Российская Федерация
2 Московский физико-технический институт (государственный университет), Долгопрудный, Московская обл., Российская Федерация
Аннотация
Испарение с поверхности почвы является важной и зачастую трудно определяемой частью модели вертикального передвижения влаги в почве. В настоящей работе задача нахождения испарения сформулирована как задача оптимального управления. Целевая функция — среднеквад-ратическое отклонение значений влажности почвы на различных глубинах, полученных в результате расчетов в рамках используемой модели, от некоторых предписанных значений. Поиск численного решения выполнен методом наискорейшего спуска, при этом градиент целевой функции рассчитан по формулам быстрого автоматического дифференцирования. Кроме того, метод быстрого автоматического дифференцирования был применен для оценки чувствительности влажности почвы на различной глубине к изменениям испарения. Полученные оценки позволили определить эффективную (небольшую) глубину подповерхностного почвенного слоя, в котором и следует рассчитывать целевую функцию. Численные расчеты показали, что это существенно ускоряет сходимость процесса поиска решения и уменьшает время его проведения
Ключевые слова
Модель вертикального передвижения влаги в почве, целевая функция, метод наискорейшего спуска, метод быстрого автоматического дифференцирования
Поступила в редакцию 04.05.2016 © МГТУ им. Н.Э. Баумана, 2016
Работа выполнена при поддержке Российского фонда фундаментальных исследований № 15-07-08952
Введение. При изучении процессов формирования стока на водосборе существенная роль отводится моделированию вертикального передвижения влаги в почве. Такая модель содержит несколько параметров, определение которых представляет собой отдельную и часто трудновыполнимую задачу. Значения большинства параметров физико-математических моделей задают либо из физических соображений, либо находят с помощью эмпирически установленных соотношений по измеряемым характеристикам речного водосбора. Поиск таких
соотношений — отнюдь не простая задача и не может быть решена для всех параметров (часть из них часто подбирают по измерениям стока). Некоторые параметры определяют экспериментально, что также бывает трудноосуществимо.
Из всевозможных способов определения параметров модели также следует упомянуть способ, когда искомые параметры находят по некоторым доступным экспериментальным данным. В качестве таких доступных данных можно рассматривать измеренные профили влажности почвы, т. е. распределение влаги в почве по глубине. Утверждение о доступности указанных данных подкрепляется следующим фактом: на европейской части Российской Федерации находятся около 500 станций, на которых на постоянной основе в мониторинговом режиме осуществляют измерения профилей влажности почвы.
Часть рассматриваемой модели — интенсивность испарения с поверхности почвы, представляющая собой линейный поток влаги через поверхность почвы в атмосферу, выраженный в единицах слоя воды в единицу времени. В научной литературе гидрометеорологического профиля интенсивность испарения часто называют испарением [1]. Испарение является одной из важнейших составляющих водного баланса речных бассейнов и отдельных участков земной поверхности. Данные об испарении очень важны для определения влагообмена между почвой и атмосферой и используются в климатологии и гидрологии. При этом испарение — одна из самых трудно определяемых величин, входящих в модель. Так, измерить испарение с поверхности почвы намного труднее, чем испарение с водной поверхности. Соответствующие приборы — почвенные испарители — существуют, но определяемое ими испарение из вырезанных монолитов почвы может отличаться от испарения в естественной обстановке. Лабораторное измерение возможно на специально оборудованных станциях, число которых на европейской части Российской Федерации не превышает 20, и проводится нерегулярно. В настоящее время не существует точных и универсальных методов расчета испарения. Большинство геофизических методов позволяют оценить испарение лишь за большие промежутки времени (недели и месяцы) и на больших территориях. Метод Пенмана позволяет вычислять испарение, исходя из сложной зависимости его от радиационного баланса, скорости ветра, температуры и влажности воздуха. Однако в этом методе используется большое количество метеорологических данных, что ограничивает его применение, так как далеко не везде на метеорологических станциях проводят все необходимые для расчетов наблюдения.
В настоящей работе рассмотрена задача нахождения испарения с поверхности почвы по профилям влажности почвы. Задача определения испарения сформулирована как задача оптимального управления, в которой управлением является искомое испарение, а целевой функцией — среднеквадратическое отклонение значений влажности почвы, рассчитанных в соответствии с моделью передвижения влаги в почве при выбранном испарении, от некоторых предписанных значений, назовем их «экспериментальными данными». Сравнение значений влажности почвы происходит в некоторой заданной области. С помощью
метода быстрого автоматического дифференцирования (БАД) оценивается чувствительность влажности на различной глубине к небольшим изменениям испарения. Эти оценки позволили определить область, в которой проводится сравнение значений влажности, вычисленных в рамках заданной модели, с предписанными значениями так, чтобы повысить эффективность оптимизационного процесса.
Постановка задачи. Предположим, что почва представляет собой изотермическую недеформируемую однородную пористую среду. При соблюдении этих предположений вертикальное передвижение влаги в почве хорошо описывается одномерным нелинейным уравнением с частными производными второго порядка параболического типа.
Рассмотрим следующую начально-краевую задачу:
^ = ^В (9)^-К(0)1, (г,г) е О;
дг дг ^ 4 ' дг
0(г,0) = ф(г), г е(0,1);
0(1,0 = у(г), г е(0,Г); (1)
-|d (е)|-К (9)
= R(t)-E(t), t е(0,T);
z = 0
emin <e(o, t )<emax, t e (o, t),
где z — пространственная координата по оси, направленной вниз от поверхности почвы; t — время; 0 — искомая влажность в точке (z, t), так называемая объемная влажность почвы, выражаемая в единицах объема воды в единичном объеме почвы (безразмерная величина); Q = ( 0, L )х( 0, Т ); ф( z), t) — заданные функции; D (e), K (e) — коэффициент диффузии и гидравлическая проводимость; emin = er +s, e max = es —s, s ^ er; R(t), E(t) — интенсивности осадков и испарения, линейные потоки влаги; 0 < E (t)< G, t e(0, Т), G — некоторая константа, G > 0.
Входящие в уравнение коэффициент диффузии D (e) и гидравлическую проводимость K (e) вычисляют по формулам Ван Генухтена [2]:
К (9) = К öS °,5 [i-(1 - S1/m )"
D (9) = Ко —1 m , s°,5-1/m Г(1 -s1/m )-m + (1 -s1/m)m -2'
V ' am (9S -9r) LV ; V ;
(2)
где К0, а, т, 0Г, 05 — некоторые параметры; 5 = (0-0г )/(05-9Г).
Назовем задачу (1) прямой задачей. Задачу нахождения интенсивности испарения Е (г), г е( 0, Г), сформулируем следующим образом. Пусть на некотором
множестве О0 с О задана функция 0 (г, г), которую назовем «эксперименталь-
ными данными». Зададим цель подобрать интенсивность Е )е[0,0 ], I е(0, Т), так, чтобы соответствующее решение прямой задачи (1) было как можно ближе к функции 0 (^, I) на множестве 0о с Q. Другими словами, задача заключается в том, чтобы найти оптимальное управление Еорг (t )е[ 0,0 ], t е( 0, Т ), и соответствующее решение 0ор1 (^, t) прямой задачи, при которых функционал
1 А 2
] = — | (Э0р{ — 0) dzdt достигал бы минимума.
2 00
Перейдем к дискретному аналогу задачи (1). Разобьем интервалы (0, Ь) и (0, Т) на I и N равных подынтервалов с концевыми точками Zi = Ш,0< I < I, и ^ = тп, 0 <п <N, где т = Т^ и к = Ь[1. Аппроксимируем прямую задачу (1) с помощью конечно-разностной схемы:
01+1 _ ßn
en+1 п; _ 0,
n+1
= pn+1 ' + 1 "i__K"n+1 _ Qn+1 _"i_1 + Kn+1
01+1 _ 01+1
+ 1/2 h i+1/2 г_1/2
0 < i < I, 0 < n < N,
h
_V2
00 =фi, 0 < i < I, 0п = \п, 1 < п < N,
где 0п, Пп+у2, Кп_у2 — значения функций 0(z,t), О(0(z,t)), К(0(z,t)) в точках (zi, 1п), ((г + 1/2)к,тп), ((г-1/2)к,тп); фг, \\п — значения функций ф(z), \(t) в точках zi и . Конечно-разностная аппроксимация второго порядка точности левого краевого условия [3] приводит к выражению [1]
01+1 дя
0 _00
=hD
-»n+1
Л/2
0Г1 _9, h
n+1
f n+1 ■_ K1/2 ■
Rn+1 _ E
n+1
0 < n < N.
Здесь Яп+1, Еп+1 — значения функций Я(t) и Е(t) в точках Р+1 =т(п + 1). В результате получаем дискретный аналог прямой задачи (1):
ф; = _
г | Dh i00
9; + 1-0; _ + - 2
2_
h D/2
9min < 0; < 9max ,
х h
1 < n < N;
R; _En ) =
) =
ф; = hi Diy2 0;_1 _
1 1 /
D1+1/2 + D1_1/2,
0;+^т D1+у 2 0;+1
(3)
011+1 / х hV
K 1/2 Ki+1/2 ^
= 0, 1 < i < 1, 1 < я < N,
Фп = \п —0п = 0, 1 < п < N, 00 = фг, 0 < г < I,
где Еп е[ 0,0 ], 1 < п < N. При этом коэффициент диффузии О и гидравлическую проводимость К в промежуточных точках будем вычислять по формулам
х
х
впвп КпКп
В" = 2 ' '+1 ; К"+1/2 = 2 К ¡+1 , 1 < п < N, 0 < г <1. (4)
¡+1/2 Вп + Вгп+1 '+1'2 Кп + Кп+1
Пусть Q0 = |(г, I): г = ¡Н, t = Iх, I е А, I е В}, где А, В — некоторые множества
такие, что А с А0, В с В0, А0 = (0,1,2,...,I}, В0 = {1,2,...,N1. Зададим целевой функционал вида
W
(0.u) = 2 "в") (5)
2 jeA ие5
Дискретная задача оптимального управления формулируется следующим образом: найти оптимальное управление иор = {ЕпР1, п = 1,.,N}, Едо е[0,С],
п = 1,...,N, и соответствующее оптимальное решение задачи (3) такие, что функционал (5) достигал бы минимального значения.
В настоящей работе исследован вопрос о выборе области Q0, в которой проводится сравнение значений влажности почвы, полученных в результате решения задачи (3) при выбранных значениях испарения, и предписанных значений влажности в соответствующих точках. Для этого оценивают чувствительность влажности почвы в точках, находящихся на различной глубине, к незначительным изменениям значений испарения. Оценивание выполняют с помощью метода БАД. Изложим кратко суть этого метода.
Необходимость вычисления точного градиента функции, значение которой находят при выполнении некоторого алгоритма, привела к созданию метода БАД [5, 6, 9]. Этот метод позволяет находить точные значения производных сложных функций, переменные которых связаны функциональными связями. В Российской Федерации метод БАД развивался в процессе разработки и совершенствования методов решения конечномерных задач оптимизации, получаемых в результате дискретизации задач оптимального управления [4, 7, 8]. Формулы БАД исследователи получали различными путями. Приведем наиболее ясный и в то же время общий способ [7] получения формул БАД для определения производных сложной функции, основанный на теореме о неявной функции.
Предположим, что для векторов г е Я" и и е Яг дифференцируемые функции Ш(г,и) и Ф(г,и) определяют отображения Ш: К" х№ ^&1 и Ф: Яп х № ^ Яп. Пусть г и и удовлетворяют системе из п скалярных алгебраических уравнений:
Ф(г, и ) = 0п, (6)
где 0п — нулевой п-мерный вектор. Предположим, что матрица Ф^ (г, и) неособенная. Тогда по теореме о неявной функции система связей (6) определяет непрерывно-дифференцируемую функцию г = г(и). Согласно методу БАД, при этих предположениях градиент функции Ш (г (и), и) вычисляют по формуле
т((М),М) = Wu ((и),и)Фт„ ((и),и)). (7)
du
Входящий в формулу (7) вектор p е Яп является множителем Лагранжа и определяется в результате решения системы линейных уравнений:
Wz (z (u), u) + Фт (z (u), u) p = 0„. (8)
Линейная система (8) является сопряженной к исходной системе связей (6).
Метод БАД был применен для вычисления производной 0п, i е А0, п е В0, по Ek, к е B0, с целью оценить чувствительность влажности 0гп к малым изменениям испарения Ek. Кроме того, поиск численного решения дискретной задачи оптимального управления проводился методом наискорейшего спуска, при этом градиент целевой функции W (0, u) рассчитывался по формулам БАД. Соотношения для вычисления градиента целевой функции W (0, u ) имеют вид
dW(0(u),u) , , . ч / / ч ч -;= Wu (0(u),u)ФU (0(u),u)p; (9)
du
We(0(u), u) + Ф0(0(u), u) p = 0М, (10)
где M = ( I + 1) N; Фт =[ф0, ф1,... , Ф1, Ф0, ф2,... , ф2,. , ФN, ФN ,---, ФN]
0т =[00,01,.,01,00,02,.,02,---,eN,eN,---,eN] (уравнениями связи являются соотношения (3)). В случае оценивания чувствительности влажности к изменению испарения полагаем в формулах (9), (10) W(0, u) = 0П , i е А0, п е B0, а в случае вычисления градиента целевой функции в указанные формулы подставляем функцию W (0, u) из (5).
Формулы (9), (10) для вычисления градиента функции W(0, ы) являются точными. При этом значение градиента находится по указанным формулам в результате некоторого вычислительного процесса. Возникает вопрос об устойчивости этой вычислительной процедуры. Формула (9) содержит операции сложения и умножения, производные функций W(0, u) и Ф(0, u), а также множитель Лагранжа p, который определяют в результате решения линейной системы уравнений (10). Учитывая вид соотношений (3), (4), можно сделать вывод о том, что система уравнений (10) расщепляется на N подсистем уравнений, каждая из которых соответствует некоторому временному слою. Двигаясь от №го слоя к первому, можно решать последовательно каждую такую подсистему уравнений отдельно от других. Легко заметить, что основные матрицы этих подсистем уравнений имеют трехдиагональную структуру:
А (') х ^—1) — С (¿)х (') + В (¿) х ^ +1) = —Р (¿), i = 0,...,! (11)
Естественно, каждая такая подсистема решается методом прогонки. Дифференцируя Ф(0, и) из (3), получаем следующие формулы для нахождения ненулевых элементов основной матрицы подсистемы, соответствующей и-му временному слою:
а (о ) = о, а (1) = А ои/2 + А (2 ) (в?-в8)-к (К?2 );
h 12 h2\ Y2/e^ 1 h
A() = FD42 + Ы ) ( -e?-1)-h^ ),
) = Dy 2 + ( 2 L (-en Ь h1 ( 2 L ;
в
(i) = D." У2 + ±Ky2 L ( -0« ) + h)e
B (i)= h2 » V2 h
В (I -1) = 0, B (I) = 0;
i = 2,..., I;
i = 1,..., I - 2;
C (0 ) = C (i ) =
1+A D
X + h2 DV2
h
-(DV2)en (e0-en) + f(2}
X + (Dr+V2 + D-1/2 )] + ^(/2) (-en-i) + (Dr+V2) ( -e+ +
(iCn+V2 L -(k-i/2 )e
1
+ —
h
, i = 1,., I -1; C (I) = 1.
С учетом приведенных выше соотношений имеем
C (0 ) = 1 + A (0 ) + 2В (0 ); х
C (1)=1 + 0,5A (1) + В (1); х
C(i)=1 + A(i) + В (i), i = 2,.,I-2;
(12)
C (I -1)=r A (I -1)+B (I -1)+F d42 + Ь (Dn-12 ) ( -1-en)+1(KKn-V2)
Функции О (в), К (в) и в/, 1 < / < N, о < I < I, ограничены в рассматриваемой области. Из формул для А (г), В (г), 0 < г < I, и соотношений (12) ясно следующее: можно подобрать такое соотношение между тик, что |С (г )|> |А (г )| + |В (г )|, т. е. все строки с нулевой по (I - 1)-ю будут обладать свойством строгого диагонального преобладания. В связи с этим вычисление всех прогоночных коэффициентов будет корректным и устойчивым.
В начале проведения обратной прогонки при вычислении х1 по формулам прогонки происходит деление на С (I), С (I ) = 1, так как В (I — 1) = 0. Поэтому и здесь не возникает осложнений несмотря на то, что !-я строка матрицы рассматриваемой системы не обладает свойством диагонального преобладания.
Таким образом, процесс вычисления множителей Лагранжа р устойчив, и, следовательно, процедура вычисления градиента по формулам БАД также устойчива.
Численные результаты. Задача рассмотрена при следующих значениях входных параметров:
Ь = 90 см, Т = 122 сут., 9тщ = 0,065 см3/см3, 9тах = 0,51 см3/см3, в = 10"10, ф(г) = 0,3, ге(0,Ь), К0 = 103,68 см/сут., а = 0,075 1/см, т = 0,471.
Зависимости влажности почвы t) при г = 90 и слоя выпавших осадков в сантиметрах за сутки Я (t) от времени приведены на рис. 1, а и б.
щ
0,200
0,150
0,010
0,005
О 20 40 60 80 100 120 t, с а
20 40 60 80 100 120 t, с б
0 20 40 60 80 100 120 г, с в
Рис. 1. Зависимости влажности почвы ) при г = 90 (а), слоя выпавших осадков в сантиметрах Я() (б), среднесуточной интенсивности испарения ЕГие () (в), соответствующей ЕГие (), от времени
Как правило, на практике наблюдаются суммарные за сутки значения осадков и испарения. Однако шаг по времени 1 сут. оказался неприемлемым в смысле точности аппроксимации непрерывного решения. Шаги по времени и по пространству дробились до тех пор, пока последующее их уменьшение не привело к тому, что вновь рассчитанные значения влажности отличались от значений, вычисленных на предыдущем этапе, приблизительно на 1-10"4. В результате шаг
по времени г составил 1/100 сут., а шаг по пространству Н — 1 см. Таким образом, в рассматриваемом случае I = 90, N = 12 200. Поскольку на практике доступными данными оказываются суммарные за сутки значения осадков и испарения, было принято допущение о постоянстве интенсивностей осадков и испарения в течение суток.
С применением формул БАД (9), (10) были рассчитаны производные влажности на разных глубинах по испарению. Из (3) следует, что производные 0", г е А0, по Ек равны нулю, если к > п, к е В0, п е В0. Численные расчеты показали, что на каждом временном слое п, п е В0, абсолютная величина производной
0п, г е А0, по Ек, к < п, естественно, убывает с увеличением г и уменьшением к. Пример результатов расчетов, относящихся к 11 810-му временному слою, приведен в таблице. Таблица состоит из 11 столбцов и 21 строки и содержит значения компонент градиентов влажностей в точках, находящихся на различной глубине и принадлежащих одному временному слою, а именно: 11 810-му временному слою. Каждая строка таблицы слева направо содержит значение индекса г, далее значения компонент градиента 0п (Е) с (п - 9 )-й по п-ю, где п = 11 810.
Наиболее значимые компоненты градиентов влажностей на 11 810-м слое
i gr[n - 9] gr[n - 8] gr[n - 7] gr[n - 6] gr[n - 5] gr[n - 4] gr[n - 3] gr[n - 2] gr[n -1] gr[n]
0 -0,006 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,01 -0,013 -0,024
1 -0,006 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,01 -0,012 -0,022
2 -0,006 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,01 -0,012 -0,021
3 -0,006 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,01 -0,012 -0,019
4 -0,006 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,009 -0,012 -0,018
5 -0,006 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,009 -0,012 -0,016
6 -0,006 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,009 -0,011 -0,015
7 -0,006 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,009 -0,011 -0,014
8 -0,005 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,009 -0,011 -0,013
9 -0,005 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,009 -0,01 -0,012
10 -0,005 -0,006 -0,006 -0,006 -0,006 -0,007 -0,008 -0,009 -0,01 -0,011
11 -0,005 -0,006 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,009 -0,01
12 -0,005 -0,005 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,009 -0,009
13 -0,005 -0,005 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,009 -0,008
14 -0,005 -0,005 -0,006 -0,006 -0,006 -0,006 -0,007 -0,008 -0,008 -0,008
15 -0,005 -0,005 -0,005 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,007
16 -0,005 -0,005 -0,005 -0,006 -0,006 -0,006 -0,007 -0,007 -0,008 -0,007
17 -0,005 -0,005 -0,005 -0,005 -0,006 -0,006 -0,006 -0,007 -0,007 -0,006
18 -0,005 -0,005 -0,005 -0,005 -0,006 -0,006 -0,006 -0,007 -0,007 -0,006
19 -0,005 -0,005 -0,005 -0,005 -0,006 -0,006 -0,006 -0,006 -0,006 -0,005
20 -0,005 -0,005 -0,005 -0,005 -0,005 -0,006 -0,006 -0,006 -0,006 -0,005
Анализ полученных результатов показал, что наиболее чувствительными к изменению испарения оказались влажности в точках, принадлежащих подповерхностному слою почвы толщиной около 2о см. В точках, расположенных на глубине 2о см, абсолютная величина производной влажности по испарению на том же временном слое иногда оказывалась меньше соответствующей величины на нулевой глубине почти в 5 раз. Как правило, абсолютная величина производной влажности по испарению на том же временном слое не превышала значения 0,005. На глубине более либо 50 см, либо 70 см влажность оказывается нечувствительной к изменению испарения. В силу этого сравнивать вычисленные значения влажности с предписанными значениями в таких точках не имеет смысла. Таким образом, чтобы сделать целевую функцию более чувствительной к изменениям испарения, следует проводить сравнение рассчитанных в рамках модели значений влажности с предписанными значениями в точках, принадлежащих области, которая находится в подповерхностном почвенном слое толщиной 20 см.
Далее расчеты проводили по следующему сценарию. На первом этапе решалась прямая конечно-разностная задача (3) с некоторой интенсивностью испарения Егие (^). График среднесуточной интенсивности испарения, соответствующей ЕГие ^), приведен на рис. 1, в.
Система (3) расщепляется на N подсистем, каждая из которых относится к г-му, г = 1,..., N, временному слою и содержит в качестве переменных значения влажности только на этом временном слое. Двигаясь в направлении от первого временного слоя к ^му временному слою, последовательно для каждого временного слоя решаем соответствующую систему уравнений отдельно от других.
На каждом временном слое коэффициенты при неизвестных зависят от этих неизвестных (см. (3)). В связи с этим решение каждой такой системы находят в результате выполнения следующего итерационного процесса. На каждой итерации значения влажности на временном слое определяют в результате решения соответствующей системы, в которой значения коэффициента диффузии и гидравлической проводимости вычислены по формулам (2), (4) с использованием значений влажности, полученных на предыдущей итерации. При известных значениях коэффициента диффузии и гидравлической проводимости система является линейной, а основная матрица системы — трехдиагональной. Система решается методом прогонки. На первой итерации значения коэффициента диффузии и гидравлической проводимости рассчитываются с применением значений влажности на предыдущем временном слое. Итерационный процесс продолжался до тех пор, пока среднеквадратическое отклонение значений влажности, полученных на предыдущей итерации, от значений влажности в текущей итерации на всем временном слое не станет менее 1 -10—7. При выбранных параметрах задачи для нахождения решения требовалось, как правило, не более пяти итераций. Полученное решение в(^,I), (2,1)е Q0, задачи (3) называлось «экспериментальными данными».
На втором этапе решали дискретную задачу оптимального управления, в которой целевая функция вычисляется по формуле (5) при А = А0 и В = В0.
Численная оптимизация проведена методом наискорейшего спуска, при этом градиент целевой функции (5) рассчитали по формулам БАД (9), (10). Шаг вдоль выбранного направления определяли в результате одномерной оптимизации функции, полученной интерполяцией целевой функции с помощью сплайнов, которые построены по 40 точкам. Итерационный процесс продолжался до тех пор, пока чебышевская норма градиента целевой функции (5) не становилась менее 1-10-14. В качестве начального управления было выбрано испарение Еш (£) = 1,0-10-2, ^е(0,Т). Задача численной оптимизации на всем промежутке времени (0, Т) допускает декомпозицию на совокупность отдельных задач оптимизации в соответствии с выбранным разбиением всего временного промежутка. В силу того, что специалисты работают с суточными данными осадков и испарения, следовало бы рассматривать множество, состоящее из 122 задач, соответствующих каждым суткам. Численные расчеты показали: для нахождения оптимального управления с требуемой точностью необходимо, как правило, не более 1000 итераций.
На третьем этапе численные расчеты были проведены по сценарию, аналогичному сценарию второго этапа, с единственным различием: целевая функция определена по формуле (5), в которой В = В0 и А = |0, 1, 2, 3}, т. е. сравнение значений влажности почвы 0( г, I), полученных в результате решения прямой задачи (3) при выбранном управлении Еп, п е В0, с предписанными значениями 0 (г, I) происходит в узлах сетки, принадлежащих подповерхностному слою толщиной 3 см. В этом случае оптимизационный процесс происходил гораздо более эффективно: во-первых, время, затрачиваемое на проведение одной итерации, заметно сократилось; во-вторых, ощутимо увеличилась скорость сходимости. Как правило, оптимизационный процесс позволял получить решение за 100-105 итераций. Найденное оптимальное управление Е0р{ отличалось от истинного управления
незначительно: \ЕШе (гп ) - Ео?х () | <
< 10,5-10-12, п = 1,...,N.
Зависимости начального, истинного управлений и найденного с помощью численных расчетов оптимального управления
Рис. 2. Зависимости начального (1), ис- от времени представлены на рис. 2. Зави-тинного (2) и оптимального (3) управ- симости истинного и найденного в резуль-лений от времени тате численных расчетов оптимального
управлений практически совпадают.
Заключение. Анализ полученных результатов показал, что идея применения метода БАД при решении задачи нахождения испарения по профилям
влажности оказалась перспективной. Во-первых, вычисление точного градиента целевой функции по формулам БАД при поиске численного решения задачи методом наискорейшего спуска позволяет восстанавливать испарение с высокой точностью. Во-вторых, применение метода БАД позволяет быстро и точно вычислять во всех точках сетки значения производных влажности по испарению. Это дает возможность выделить эффективный подповерхностный почвенный слой, в точках которого следует проводить сравнение вычисленных значений влажности с предписанными значениями, что позволяет уменьшить число точек, в которых рассчитывают целевую функцию, и повысить ее чувствительность к изменениям испарения. Выделение такого эффективного множества точек приводит вследствие уменьшения их числа к существенному уменьшению времени, затрачиваемому на проведение одной итерации, и, кроме того, повышает чувствительность целевой функции к изменениям испарения. Перечисленные факторы повышают эффективность оптимизационного процесса, увеличивая скорость сходимости к решению и уменьшая время, затрачиваемое на нахождение решения.
ЛИТЕРАТУРА
1. Кучмент Л.С., Демидов В.Н., Мотовилов Ю.Г. Формирование речного стока. М.: Наука, 1983.
2. Van Genuchten M.Th. A closed form equation for predicting the hydraulic conductivity of unsaturated soils // Soil. Sci. Soc. Am. J. 1980. Vol. 44. P. 892-898.
3. Самарский А.А. Введение в теорию разностных схем. М.: Наука, 1971.
4. Айда-Заде К.Р., Евтушенко Ю.Г. Быстрое автоматическое дифференцирование на ЭВМ // Математическое моделирование. 1989. Т. 1. № 1. С. 120-131.
5. Griewank A. On automatic differentiation // Mathematical Programming: Recent Developments and Applications. Ed. by M. Iri, K. Tanabe. Tokyo: Kluwer Academic Publishers, 1989. P. 83-108.
6. Automatic differentiation of algorithms. Theory, implementation and application. Ed. by A. Griewank, G.F. Corliss. Philadelphia: SIAM, 1991.
7. Evtushenko Yu. Automatic differentiation viewed from optimal control theory // Automatic differentiation of algorithms. Theory, implementation and application. Ed. by A. Griewank, G.F. Corliss, Philadelphia: SIAM, 1991. P. 25-30.
8. Evtushenko Yu. Computation of exact gradients in distributed dynamic systems for optimal control problem // Optimization Methods and Software. 1998. Vol. 9. P. 45-75.
9. Griewank A. Evaluating derivatives. Philadelphia: SIAM, 2000.
Дикусар Василий Васильевич — д-р физ.-мат. наук, профессор, ведущий научный сотрудник Вычислительного центра им. А.А. Дородницына ФИЦ «Информатика и управление» РАН (Российская Федерация, Москва, 119333, ул. Вавилова, д. 40).
Засухин Сергей Владимирович — аспирант Московского физико-технического института (государственного университета) (Российская Федерация, 141700, Московская обл., Долгопрудный, Институтский переулок, д. 9).
Просьба ссылаться на эту статью следующим образом:
Дикусар В.В., Засухин С.В. Применение быстрого автоматического дифференцирования при нахождении испарения с поверхности почвы // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2016. № 6. С. 42-55. DOI: 10.18698/1812-3368-2016-6-42-55
APPLICATION OF FAST AUTOMATIC DIFFERENTIATION IN THE SOIL SURFACE EVAPORATION PROBLEM
V.V. Dikusar1 S.V. Zasukhin2
dikussar@yandex.ru s.zasukhin@yandex.ru
1 Dorodnitsyn Computing Centre, Federal Research Centre Computer Science and Control, Russian Academy of Sciences, Moscow, Russian Federation
2 Moscow Institute of Physics and Technology (State University), Dolgoprudny, Moscow Region, Russian Federation
Abstract
Evaporation from the soil surface is an important and often hard-determined part of modeling the vertical water transfer in soil. The study considers the problem of determining evaporation and formulates it as an optimal control problem. The objective function is mean-square deviation of modeled values of soil moisture at various depths from some prescribed values. We found the numerical solution by the steepest descent method, and we calculated the gradient of the objective function by formulas of fast automatic differentiation (FAD). Moreover, we applied FAD method to estimate the sensitivity of soil moisture at various depths to changes of evaporation. The findings of the research made it possible to determine the effective subsurface soil layer in which it is advisable to compute the objective function. Numerical results showed that this subsurface soil layer determination has accelerated the convergence of finding solutions process and has reduced its run time
Keywords
Vertical water transfer in soil, objective function, steepest descent method, fast automatic differentiation method
REFERENCES
[1] Kuchment L.S., Demidov V.N., Motovilov Yu.G. Formirovanie rechnogo stoka [River runoff formation]. Moscow, Nauka Publ., 1983.
[2] Van Genuchten M.Th. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil. Sci. Soc. Am. J., 1980, vol. 44, pp. 892-898.
[3] Samarskiy A.A. Vvedenie v teoriyu raznostnykh skhem [Introduction to the theory of difference schemes]. Moscow, Nayka Publ., 1971.
[4] Ayda-Zade K.R., Evtushenko Yu.G. Fast automatic differentiation on computers. Matematicheskoe modelirovanie [Matem. Mod.], 1989, vol. 1, pp. 121-139.
[5] Griewank A. On automatic differentiation. Mathematical Programming: Recent Developments and Applications. Ed. by M. Iri, K. Tanabe. Tokyo, Kluwer Academic Publishers, 1989, pp. 83-108.
[6] Automatic differentiation of algorithms. Theory, implementation and application. Ed. by A. Griewank, G.F. Corliss. Philadelphia, SIAM, 1991.
[7] Evtushenko Yu. Automatic differentiation viewed from optimal control theory. Automatic differentiation of algorithms. Theory, implementation and application. Ed. by A. Griewank, G.F. Corliss, Philadelphia, SIAM, 1991, pp. 25-30.
[8] Evtushenko Yu. Computation of exact gradients in distributed dynamic systems for optimal control problem. Optimization Methods and Software, 1998, vol. 9, pp. 45-75.
[9] Griewank A. Evaluating derivatives. Philadelphia, SIAM, 2000.
Dikusar V.V. — Dr. Sci. (Phys.-Math.), Professor, leading researcher of the Dorodnitsyn Computing Centre, Federal Research Centre Computer Science and Control, Russian Academy of Sciences (ul. Vavilova 40, Moscow, 119333 Russian Federation).
Zasukhin S.V. — post-graduate student of the Moscow Institute of Physics and Technology (State University) (Institutskiy pereulok 9, Dolgoprudny, Moscow Region, 141700 Russian Federation).
Please cite this article in English as:
Dikusar V.V., Zasukhin S.V. Application of Fast Automatic Differentiation in the Soil Surface Evaporation Problem. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Estestv. Nauki [Herald of the Bauman Moscow State Tech. Univ., Nat. Sci.], 2016, no. 6, pp. 42-55. DOI: 10.18698/1812-3368-2016-6-42-55
