CC BY
52
УДК 519.677:519.688
Поиск глобального минимума в задачах параметрической идентификации*
Л.А. Хворова, К.А. Немчикова, Д.П. Ломиворотов Алтайский государственный университет (Барнаул, Россия)
Searching for a Global Minimum in Parametrical Identification Problems
L.A. Khvorova, K.A. Nemchikova, D.P. Lomivorotov Altai State University (Barnaul, Russia)
Цель исследования — разработка и апробация алгоритмов параметрической идентификации модели продукционного процесса сельскохозяйственных культур Agrotool к условиям Алтайского края, определение допустимых границ изменения параметров модели, оценка применимости модели для прогноза урожайности зерновых культур.
Рассматривается задача адаптации модели продуктивности Agrotool к условиям Алтайского края. Решается проблема идентификации параметров модели по агрометеорологическим данным региона, осуществляется исследование чувствительности модели к вариациям и точности задания входящих в нее параметров.
Для нахождения оптимальных значений параметров применяются методы глобальной оптимизации, а именно диагональный подход к решению многомерных задач. Приводится схема диагонального алгоритма. Анализируются результаты идентификации блоков водного режима, фенологического развития и продуктивности посева. Сформулированы критерии точности задания области допустимых значений параметров модели.
Ключевые слова: модель, продукционный процесс, идентификация параметров, глобальный минимум, методы глобальной оптимизации, диагональный подход, диагональный алгоритм.
DOI 10.14258/izvasu(2014)1.2-22
The goal of the research is to design and to evaluate parameter identification algorithms for Agrotool productivity model adapted for the Altai region. It is required to determine permissible limits of the model parameters changes and to perform the model validation for predicting crop yield.
In the paper, problems of Agrotool productivity model adaptation for conditions of the Altai region are discussed. Namely, there are identification of model parameters in accordance with agro-meteorological data of the region and study of the model sensitivity to variations and accuracy of model parameters.
Optimal values of parameters are evaluated by methods of global optimization with diagonal approach to a solution of multidimensional problems. The scheme of the diagonal algorithm is provided. Results of identification for blocks of soil moisture regime, phenological development, and crops productivity are analyzed. Also, the accuracy criteria for setting the permissible limits of the model parameters are elaborated.
Key words: model, production process, identification of parameters, global minimum, diagonal approach, diagonal algorithm.
1. Постановка задачи параметрической идентификации. Пусть 3= 3 (X, Р, L) — упрощенный образ системы (модель), хх е X , х = 1,пх , — совокупность входных переменных; s¡ е 5 , х = 1, п5 — совокупность переменных состояния модели; р{ е Р, х = 1,пр — совокупность параметров модели;
еЕ , х = 1, па — совокупность внутренних связей в модели между переменными (структура модели). Функция Ь = {ь1,..., Ьщ} — разрешающий оператор совокупности математических соотношений, позволяющий по заданным входам хх е X, х = 1, пх находить функции е 5 , х = 1, п5 на интервале t0 < t < :
* Работа выполнена при финансовой поддержке благотворительного фонда В. Потанина.
S(t +1) = L(X(t),S(t), P,). Данная зависимость -закон функционирования модельной системы 3.
Задача параметрической идентификации сводится к оцениванию параметров p е P, i = 1,np . Решение поставленной многомерной задачи достигается методами глобальной оптимизации [1] и заключается в следующем:
Z(P ) = min S(P)-s
PeD * r
(1)
где Р* — вектор оптимальных значений параметров, 5(Р) — переменные состояния модели, $геа1 — фактические значения переменных состояния,
В = [а, Ь] = {Р е Ея : а(]) < Р() < Ь(Д 1 < ] < пр} . (2)
Будем предполагать, что целевая функция (1) является многоэкстремальной, недифференцируемой, заданной в форме черного ящика и удовлетворяющей в области поиска В С Мп условию Липшица с неизвестной константой Липшица 0 <Ь :
|z(p 7) - z(p 7/)| < l||p '-p '
P' ,P" eD, 0<L
(3)
где а, Ь е Мп — заданные векторы, • — евклидова норма.
Определение оптимальных значений параметров является задачей поиска глобального минимума функционала (1) — невязки (расхождение между экспериментальными данными и расчетными).
Поиск оптимальных значений параметров модели представляет собой достаточно сложную и трудоемкую задачу. Определение оптимальной комбинации параметров аналитически не представляется возможным в силу сложности модельной структуры
Agrotool: блочного характера описания процессов, взаимосвязи и взаимозависимости информационных потоков внутри модели, поэтому возникает необходимость применения численных методов.
В свою очередь, проблема численного решения задачи оптимизации (1)-(3) сопряжена со значительной размерностью вектора идентифицируемых параметров модели Р, многоэкстремальностью и недифференцируемостью целевой функции (1). Необходимо также отметить, что проведение испытания целевой функции требует значительных затрат.
На рисунке 1 приведена двумерная многоэкстремальная целевая функция невязки для влажности почвы. Поведение целевой функции отличается очень узкими и глубокими пиками и впадинами, которые могут находиться между точками вычислений. Для успешного решения оптимизационной задачи (1)-(3) использование методов нелинейной локальной оптимизации оказывается недостаточным в силу наличия нескольких локальных минимумов, имеющих разные значения целевой функции. «Локальные методы, как правило, оказываются не в состоянии покинуть зоны притяжения локальных оптимумов и, следовательно, упускают глобальный экстремум» [1, с. 11].
2. Методы и алгоритмы глобальной оптимизации. Если минимизируемая функция имеет больше одного экстремума и нужно найти наименьший из них, то ни один из классических методов и алгоритмов поиска локального экстремума не решит поставленную задачу. Ввиду высокой сложности таких задач и для быстрого поиска глобального решения многоэкстремальной задачи используют эффективные алгоритмы глобального поиска, например, диагональный подход [1] или метод имитации отжига [2, с. 85-105].
При использовании диагонального подхода для целей идентификации в (2) были определены три
Рис. 1. Двумерная многоэкстремальная целевая функция невязки для влажности почвы
гиперинтервала Di (для параметров трех блоков модели: динамики влажности почвы, фенологического развития и продуктивности посева). Каждый D¡ разбивается на множество гиперинтервалов Dki, х = 1,2,3 — для параметров трех блоков модели, 1 < к < Ы{ — число гиперинтервалов в каждой области D¡. Целевые функции ^ (Р), х = 1,2,3, вычисляются в двух вершинах, а{ и Ь, главной диагонали каждого гиперинтервала Dki. В работе [1] описаны различные диагональные стратегии разбиения и приведены вычислительные схемы диагональных алгоритмов, которые авторы статьи применили в своем исследовании.
3. Результаты параметрической идентификации с помощью оптимизационных процедур. В результате проведенного исследования построена оптимизационная процедура параметрической идентификации блоков модели: динамики влажности почвы, фенологического развития и продуктивности посева [3, с. 470-472; 4, с. 171-175].
В каждой области D¡, х = 1,2,3, была построена неравномерная сетка, уплотняющаяся в окрестности глобального минимума. Так, координаты точки глобального минимума при идентификации параметров блока динамики влажности в почве определены при помощи вычисления значений невязок на сетках с различными шагами для К/
и = 100, к„ = 10, „ , к„ = 0,1 и С :
К ' К/ ' К/ ' К/
кС = 0,1, кС = 0,05 .
полнения 1920 вычислительных процедур, прежде чем был найден глобальный минимум, k — индекс текущего разбиения области D.
1. Результаты идентификации блока влагоперено-са в почве. В основу модели влагопереноса в модели Agrotool [5, с. 99-105] положено уравнение Ричардса
.Мхо =4 ^ (Р.) -1|-/(хД),
дt дх I 5 дх
где t — время; х — пространственная координата; в — объемная влажность почвы; Р — капиллярно-сорбционный потенциал почвенной влаги; к™ (Р.) — функция влагопроводности: к™ (Р5) = К/ ■ (—Р5 )С К/— коэффициент фильтрации (см/сут), С—эмпирический параметр; /(х,0 — функция стока.
Коэффициент фильтрации К/и показатель степени С определяются диагональным методом глобальной оптимизации. Целевая функция (1) принимает вид (4)
(К/ ,С) = ЕЕ|(в5о,г М-вга М)
¡=1 }=1
- ш1п . (4)
К/, СеР
Разбиение двумерного (суженного) гиперинтервала
D1k = {4,0 < К/ < 6,0, 1,1 < С < 1,9} потребовало вы-
Здесь вгеа1 (х,}) — фактические значения влагозапаса, вюи а,}) — расчетные значения, х = 1, т — номер года, т — общее число лет, за которые производится компьютерный эксперимент, } = 1, к — число фактических замеров влагозапаса в почве в течение т лет.
Для иллюстрации точности расчетов, достигнутых в результате оптимизационной процедуры (4), на рисунке 2 представлены экспериментальные и расчетные данные по динамике влагозапаса в почве под посевом яровой пшеницы.
Рис. 2. Динамика влагозапасов (ОПХ им. Докучаева, 2005 г.): сплошная линия — влагозапас, рассчитанный по модели;
■ — экспериментальные значения
8.00 * 7.00 л 6.00 2 5 00
О
и 4.00 о 3.00
О
2.00
О
с 1.00 0.00
* Л ¿° ¿ь
V V V V
<: <=■> <> ^ ^
Коэф.фильтрации, см/сут —♦— Сред. погрешность
Рис. 3. График средней относительной погрешности вычисления запасов влаги в почве в зависимости от величины
коэффициента фильтрации, С = 1,1 (суглинок тяжелый)
8.50 8.00
7.50
6.00
£ ^ ^ # # # <§> <<? к4 к"? & ^
<о <о
Коэф. фильтрации, см/сут
-Сред. погрешность
Рис. 4. График средней относительной погрешности вычисления запасов влаги в почве в зависимости от величины
коэффициента фильтрации, С = 1,4 (суглинок средний)
Для тяжелосуглинистых почв минимальная средняя погрешность по формуле (4) — 5,4% достигнута при оптимальных значениях К=5,9 и С=1,1. График средней относительной погрешности приведен на рисунке 3.
Для среднесуглинистых почв минимальная погрешность составила 7% при оптимальных значениях К = 61,2 и С = 1,4. График средней относительной погрешности приведен на рисунке 4.
Определение допустимых границ изменения параметров блоков осуществлялось с помощью методов теории чувствительности (метод прямого моделирования и метод оценок вариаций). В процессе исследования определены функции чувствительности модели как функции влияния изменений параметров на решение задачи [6, с. 128-132]. Показано, что для тяжелосуглинистых почв допустимый интервал изменения К — (4.0-6.0); значения С не только сильно влияют на динамику влажности почвы, но и на величину урожая и поэтому требования к величине С достаточно жесткие: С = 1.1. Аналогичная ситуация наблюдается и для легкосуглинистых почв: К— (40-65), С = 1.4. Анализ на чувствительность модели осуществлен также к другим гидрофизическим параметрам почвы и начальному состоянию модели.
2. Результаты идентификации блока фенологического развития. Задачами этого блока являются расчет
так называемого физиологического времени, измеряемого в градусо-днях, и сроков наступления фенологических фаз. В модели Agrotool приращение физиологического времени в день k вычисляется по формуле
дт( к) = дт0(1 -дт0/с,) ф
где
Дто(к) =
(Т„ (к) - То) при Тт (к) > то, 0 при Т„(к)<ТО,
Тт(к) — среднесуточная температура воздуха в день k, к = 1,365 ; Тд — биологический нуль, с — константа; ф: — потенциал воды в почве;
Б^ф) =
1 если ф > ,
1 + («о-1)
Ф, - Ф0
орг
Ф„ - ф0
, если Ф < фор(,
орг
фп — потенциал воды в почве, соответствующий влажности завядания.
Величина прироста биологического времени определяется по формуле
т(к) = ]т дт( 1),
¡=к„
7.00
г 6.50
где kg — номер дня сева; k — номер текущего дня. Очередная фаза наступает при достижении величиной т(к) некоторого порогового значения ТРк (1Рк), зависящего от порядкового номера фазы 1РЬ: т(к) > ТРЬ (1Рк).
При идентификации пороговых значений определялось минимальное расхождение между расчетными и фактическими датами наступления фенофаз. В таблице 1 приведены оптимальные значения биологических порогов.
Таблица 1
Результаты идентификации пороговых значений
Таблица 2
Расчетные и экспериментальные величины урожаев*, ц/га
Название фазы Пороговые значения
Всходы 82
Кущение 176
Выход в трубку 214,2
Колошение 363,4
Цветение 420
Молочная спелость 540
Восковая спелость 700
Полная спелость 805
Урожайность
Годы Фактическая Расчетная
1999 12,2 13,2
2000 22,2 22,0
2001 20,2 19,6
2002 22,1 23,3
2003 15,4 14,1
2004 14,1 10,2
2005 14,1 19,0
2006 25,4 20,4
2007 17,6 17,4
2008 18,9 20,1
2009 35,5 38,5
2010 19,4 21,6
3. Результаты идентификации блока продуктивности растений. Результатом окончательной идентификации параметров модели является величина урожайности культуры. Для ОПХ им. Докучаева Алтайского края расчетные и фактические величины урожаев яровой пшеницы представлены в таблице 2.
Численные эксперименты с использованием описанных оптимизационных процедур поиска глобального минимума в задаче идентификации параметров модели и анализа на чувствительность позволили: 1) разработать критерии точности задания областей допустимых значений параметров модели; 2) выявить особенности в настройке параметров блоков
* Средняя относительная погрешность составила 12%.
(среди почв с текстурой суглинок средний и суглинок тяжелый выделены почвы с текстурой суглинок легкий, что существенно сказалось на качестве поиска глобального минимума. Параметры блока фенологического развития зависят от типа почвы и сорта культуры (в модели такие зависимости не предусмотрены); параметры продуктивности посева напрямую зависят от сорта культуры (в модели параметры одинаковые для всех сортов). Учет особенностей в настройке параметров показал достаточно высокую эффективность описанного оптимизационного подхода к поиску глобальных минимумов в задаче идентификации параметров и позволил дать высокую качественную оценку применимости модели Agrotool для прогноза урожайности зерновых культур в условиях Алтайского края.
Библиографический список
1. Сергеев Я.Д., Квасов Д.Е. Диагональные методы глобальной оптимизации. — М., 2008.
2. Шарый С.П. Стохастические подходы в интервальной глобальной оптимизации // Труды XIII Байкальской Междунар. школы-семинара «Методы оптимизации и их приложения». Т. 4. Интервальный анализ. — Иркутск, 2005.
3. Хворова Л.А. Идентификация параметров модели фенологического развития зерновых культур в условиях Алтайского края // Обозрение прикладной и промышленной математики. — 2010. — Т. 17, вып. 3.
4. Хворова Л.А. Оптимизация процесса структурно-параметрической идентификации моделей продуктивности агроэкосистем // Известия Алт. гос. ун-та — 2012. — №1/1(73).
5. Хворова Л.А., Топаж А.Г. Построение моделей аг-роэкосистем и их адаптация к конкретным условиям // Научно-технические ведомости СПбГПУ — 2011. — №1(115).
6. Хворова Л.А. Методы исследования чувствительности моделей продуктивности агроэкосистем // Известия Алт. гос. ун-та — 2013. — №1/1(77).
