CC BY
55
УДК 519.6
Л. А. Хворова
Методы исследования чувствительности моделей
•к
продуктивности агроэкосистем
L. A. Khvorova
The Research on the Sensitivity of Models of Agroecosystems’ Productivity
Рассматриваются методы исследования чувствительности моделей и их приложение к динамическим имитационным моделям продукционного процесса сельскохозяйственных культур. Обсуждается практическая реализация этих методов для оценок вариаций выходных переменных и переменных состояния модели Agrotool. Ключевые слова: агроэкосистема, продукционный процесс, модель, чувствительность моделей к вариациям параметров.
Проблема решения практических задач принятия решений в агроэкологии с помощью математических моделей продуктивности агроэкосистем не исчерпывается разработкой универсальных моделей и эффективных численных методов для их компьютерной реализации. Необходимо иметь представление о качестве модели, степени ее соответствия реальной физической системе, устойчивости и чувствительности к вариациям входных параметров. Последнее является предметом исследования в данной статье и связано с тем, что значения некоторых параметров, например входных начальных данных, вычисляются по данным измерений, содержащим как случайные, так и систематические ошибки. Кроме того, при численном моделировании каждой дискретной модели присущи ошибки аппроксимации, округления, погрешности за счет итерационных процессов и др.
Содержание данной статьи составляют методы теории чувствительности для дискретных моделей продукционного процесса культурных растений. Обсуждается практическая реализация этих методов для оценок вариаций выходных переменных и переменных состояния модели Agrotool. Исследование чувствительности модели Agrotool является необходимым этапом при адаптации модели к новым условиям и культурам и положительно влияет на качество прогнозирования состояния агроэкосистемы и при выработке агротехно-логических решений в процессе формирования урожая.
Реализованная в виде программного кода динамическая модель продукционного процесса представляет собой формализованный алгоритм рекуррентного пошагового пересчета вектора состояния динамических характеристик агроэкосистемы:
The article considers methods of investigating the sensitivity of the models and their application to the dynamic simulation models of agricultural production. Practical realization of these methods for estimating variations of Agrotool model’s target variables and variables of condition is discussed.
Key words: agroecosystem, production process, model, sensitivity of the models to variations of the parameters.
S (k + 1) = L (X (k), S (k), P Z, U (k)), (1)
S (0) = S0, k = 0,1,-,T, (2)
где S (k), S (k + 1) — векторы переменных состояния системы на k-м и k+1 -м шагах; X (k) — совокупность входных переменных, в состав которых входят и неконтролируемые внешние воздействия (погода); P — вектор статических параметров модели; Z — совокупность внутренних связей в модели между переменными — структура модели; U (k) — вектор управляющих воздействий (агротехника); функция L есть не что иное, как разрешающий оператор совокупности математических соотношений, позволяющих по заданным входам (совокупности внешних воздействий) X с той или иной определенностью находить функции S ; S0 — начальное состояние системы (начальное условие); T — время окончания процесса моделирования.
Для счета модели задают значения вектора параметров P, внешние воздействия X и начальное состояние системы Sg. Вектор состояния S—это совокупность величин, исчерпывающим образом характеризующих состояние моделируемой системы в момент времени k. К переменным состояниям относятся температура и влажность почвенных слоев, биомасса отдельных органов растения и др. В качестве неконтролируемых внешних воздействий в модели выступают суточные погодные данные. Под контролируемыми воздействиями понимаются величины, характеризующие агротехнологию: сроки сева и нормы высева, сроки и нормы внесения удобрений, сроки и нормы полива и др.
1. Методы исследования чувствительности. Методы исследования чувствительности предназначены для оценки влияния вариаций входных парамет-
* Работа выполнена в рамках государственного задания «Изучение процессов конвекции и теплопереноса в анизотропных областях и областях с границами раздела» №7.3975.2011.
ров на вектор состояния и выяснения роли различных факторов при формировании моделируемых процессов. В настоящее время наиболее распространенным методом изучения чувствительности моделей является метод прямого моделирования [1], основная идея которого состоит в следующем. Задача решается при невозмущенном р = Р є D и возмущенном Р2 = Р + ЗР є D значениях вектора параметров, где Б — область допустимых значений параметров. В результате получаются два значения вектора состояния, которые, используя (1), можно предоставить следующим образом:
S1 (к + 1) = L (X (к), S1 (к), Рр 2, и (к)) и S2 (к + 1) = L (X (к), S2 (к), Р, 2, и (к)). (3)
Искомая вариация получается по формуле
ЗБ = Б1 - S2. (4)
Для проведения численных экспериментов метод прямого моделирования является наиболее простым и универсальным. Его целесообразно применять при оценке влияния больших вариаций [1].
При малых вариациях ЗР и ЗБ для нахождения ЗБ лучше использовать уравнения в вариациях. Это линейные уравнения, которые получаются в результате линеаризации исходной нелинейной системы в окрестности невозмущенных состояний и параметров. Поскольку для компьютерной реализации задача должна быть записана в дискретной форме, то и уравнение для вариаций удобно вывести из (1). Пусть Р0 — невозмущенное значение вектора параметров, а Б( Р0) — решение задачи (1) при этих параметрах. Представим вектор параметров в окрестности Р0 в виде:
Р = Р0 +%-ЗР, Р,Р0 є D, (5)
где (Ц — вещественный параметр, и предположим, что и вектор состояния можно представить в таком же виде:
Б = Б(Р0) + %-ЗБ . (6)
Тогда, используя (5), (6) и проводя формально операции дифференцирования функционала в тождестве (1), получим систему уравнений в вариациях. Уравнения в вариациях, в силу линейности полученной задачи относительно вариаций векторов состояния и параметров, при специальном задании ЗР позволяют оценивать непосредственно функции чувствительности вектора состояния к вариациям параметров. Как будет показано далее, эти функции определяются как частные производные компонент вектора состояния по компонентам вектора параметров.
Для оценок функций чувствительности применим и метод прямого моделирования. Вычисления организуются по схеме (3), (4). При каждом решении системы уравнений (1), (2) одной из компонент Р вектора Р дается небольшое приращение ЗРІ относительно невозмущенного значения. Полученная в результате
операций (3), (4) вариация вектора состояния, которую обозначим через ЗS(,), делится на ЗР1. Таким образом находится приближенное значение функции чувствительности к вариациям параметра Р .
Если требуется оценить вариации функционалов, определенных на множестве решений задачи, то в этих случаях функционал можно принять в качестве меры чувствительности модели к вариациям входных величин. Примерами таких функционалов могут быть полная энергия системы, среднее значение какого-либо компонента вектора состояния по области Б или по ее части и др. В зависимости от целей исследования функционал может быть выбран многими способами. Рассмотрим метод оценки вариаций для функционала общего вида, предполагая при этом его ограниченность, непрерывность и дифференцируемость на множествах функций 5 и Р [2]. Для определенности будем считать, что функционал явно зависит от компонент вектора состояния и неявно от компонент вектора Р . Обозначим его через 10 (5).
Из сделанных предположений относительно дифференцируемости функционала 10 (5) следует, что существуют векторы
\31о (Б)
gradSI (Б) =
дБ
(7)
gradрI (Б) =
гад
дР
Р є D
где символ частной производной обозначает операцию дифференцирования скалярной функции по векторному аргументу, результатом которой является вектор-градиент этой функции в направлении вектора-аргумента.
Определим вариации ЗІ0 (Б) функционала 10 (Б), возникающие при варьировании вектор-функций Б и Р . С учетом (7) выражение для вариации ЗІ0 (Б) при варьировании вектор-функции Б запишется следующим образом:
ЗІ0 (Б) = (gmdsI0 (Б ),ЗБ) =
дІ0 (Б) дБ
,ЗБ I.
(9)
Учитывая (8), выражение для З10 (Б) при варьировании вектор-функции Р можно записать в виде:
ЗІ0 (Б) = (gradрIo (Б),ЗР).
(10)
Скалярные произведения в (9), (10) определяются структурой функционала З10 (Б) и соотношением (1) таким образом, чтобы все операции в них имели смысл. Это замечание относится к тому, что компоненты вектор-функций 5 и Р имеют различные по физическому смыслу размерности, и при выборе функционалов и скалярных произведений необходимо учитывать это различие.
Ввиду явной зависимости функционала 10 (Б) от 5 вычисление вариации по формуле (9) производится обычным образом. Процедура вычислений опреде-
ляется конкретным видом функционала и реализует выражение
д
SI0(S) = (gradsI0(S),SS) = lim—10(S + ^-SS), (11)
дд
где д — некоторый вещественный параметр. Компоненты вектора gradSI0 (S) вычисляются непосредственно по формуле (7) как частные производные скалярной функции многих переменных по одной из них. Поскольку функционал I0 (S) в общем случае нелинейный, то удобнее определять вектор gradSI0 (S), исходя из соотношения (11), следующим образом:
gradsI0S (S) =
д
дSS
lim—10 (S + ZSS)
^0 дд °V
(12)
Уравнение (10) является основным соотношением теории чувствительности моделей. Оно связывает вариации функционалов с вариациями входных параметров. Компоненты вектора gradSI0 (Б) характеризуют вклад соответствующих им вариаций компонент вектора Р в вариацию функционала, и поэтому они называются функциями влияния вариаций параметров на функционал или функциями чувствительности модели к вариациям входных величин, если в качестве меры чувствительности взят функционал 10 (Б) . Вычисление вариаций функционала 10 (Б) по формуле (10) в зависимости от вариаций вектора Р и компонент вектора gradSI0 (Б) представляет более сложную задачу, чем вычисление по формулам (11), (12), так как при этом необходимо учитывать неявную зависимость между векторами £ и Р . Чтобы учесть эту зависимость, используют метод множителей Лагранжа, который связан с теорией сопряженных уравнений [1, 2].
2. Определение функций чувствительности. Определим функции чувствительности модели (1), (2) как функции влияния изменений параметров на решение задачи. В таком случае исследование чувствительности относится к изучению поведения модели в пространстве параметров. Математически определение функций чувствительности сводится к вычислению частных производных от искомого решения или от функционалов, определенных на множестве решений, по параметрам модели в окрестности невозмущенных значений [1]. Предполагается, что все составляющие вектора параметров Р известны.
Для целей практического использования моделей возникает вопрос о требуемой точности определения параметров. Эта задача может быть решена путем исследования чувствительности модели к вариациям той или иной составляющей вектора Р. Если какая-либо из выходных переменных или переменных состояния модели значительно изменяется даже при малой вариации параметра, то к точности определения этого параметра должны предъявляться особые требования. При незначительном же изменении выхода или состояния модели при вариации данно-
го параметра он может быть оценен приближенно. Обозначим эту варьируемую составляющую вектора Р через р. Пусть Р1 — часть вектора Р, включающая неизменяемые параметры. В таком случае уравнения модели (1), (2) перепишутся в виде:
Б ( (к + 1), р) = L (X (к), Б (к, р), Р1, р, 2, и (к)), (13)
Б (0) = Б0 (р), к = 0,1,...,Т. (14)
В (14) отражен тот факт, что в общем случае от варьируемого параметра может зависеть и начальное условие. Обычно параметры входят в уравнения нелинейно. Поэтому для исследования зависимости поведения модели от параметра вблизи его некоторого значения система (3), (4) линеаризуется. Продифференцируем обе части равенств (13), (14) по параметру р:
дБ ((k +1), p) дp
дЬ
где — =
дS
дL дS(k, p) ді дS дp дp
дS(0) = дS0 дp дp
дS,
— матрица размера
— и-мерные векторы.
1 'I,1=1,
дБ_ дL(к,р) дБ(0) др др др
Частные производные от Б (к,р) по параметру р называются функциями влияния параметрар или функциями чувствительности по отношению к вариациям параметра р. Обозначая их через V (к,р), получим:
пл,д^ дL дS0
v(k +1, p) = — • v(k, p) + — , v(0, p) = — до дp дp
(15)
Соотношения (15) называются уравнениями чувствительности. Очевидно, что они являются линейными. Для вычисления функций чувствительности необходимо иметь решение исходной системы (1), (2) при базовом значении параметра р.
Поскольку объектом исследования является компьютерная версия модели, производные в соотношениях (15) должны быть заменены конечными приращениями:
Д5 = R Ар
(16)
где 5 — произвольная характеристика системы; р — параметр из некоторой области допустимых значений Б; я0 и р0 — заданные базовые значения исследуемой переменной и параметра, которые называются невозмущенными; Ар — вариация (приращение — положительное или отрицательное) параметра; Дя — изменение исследуемой величины, вызванное вариацией параметра р; р1 = р0 + Др — возмущенное значение параметра. В соотношении (16) безразмерный коэффициент R характеризует чувствительность: чем больше окажется значение показателя R, тем более жесткие требования предъявляются к точности задания параметра р.
Если задаться определенной границей по точности расчета переменной 5 в модели, например 10%, то легко вычислить требуемую точность параметра:
ар«
0.1
р0 К
Все компьютерные эксперименты по анализу на чувствительность выполнены в имитационной системе Agrotool [3]. Поэтому необходимо отметить, что исследование чувствительности модели на самом деле связано с исследованием чувствительности алгоритма моделирования.
3. Исследование чувствительности к почвенно-гидрологическим параметрам. Почва в модели рассматривается как среда, неоднородная по своим гидрофизическим свойствам. Гидрофизические свойства каждого слоя почвы характеризуются следующими константами: максимальной гигроскопичностью (МГ), влажностью завядания (ВЗ), наименьшей влагоемкостью (НВ) и влажностью насыщения (ПВ). Именно они определяют форму зависимости потенциала почвенной влаги Р от влажности почвы 0, являющейся основной гидрофизической характеристикой (ОГХ) почвы:
в(Р) =М •
(Р) 1+Ь(-Р)а
О
к"
Ї
и
к
5
к
я
X
и
е-
с
с
55
К
я
ю
а
(17)
Показатели МГ и ПВ входят в соотношение (17) непосредственно, а две эмпирические константы а и Ь определяются по значениям ВЗ и НВ [4]. Вид зависимости потенциала почвенной влаги при различных значениях гидрофизических констант приведен на рисунке 1.
Й Объемная влажность почвы, см3/см3
Рис. 1. Влияние почвенно-гидрофизических параметров на кривую ОГХ: 1 — кривая ОГХ с базовыми значениями параметров; 2 — кривая ОГХ с измененным параметром НВ + 10%; 3 — кривая ОГХ с измененным параметром ВЗ + 10%; 4 — кривая ОГХ с измененным параметром МГ + 10%; 5 — кривая ОГХ с измененным параметром ПВ + 10%
Как видно из рисунка, при увеличении НВ кривая водоудерживания смещается вверх по сравнению с исходной при различных значениях влажности почвы. При увеличении ВЗ кривая проходит выше исходной при низких значениях влажности почвы и ниже исходной при высоких величинах влажности. При увеличении ПВ картина становится обратной, а при варьировании МГ вид кривой фактически не изменяется.
е;
(С
Дни от момента сева
Рис. 2. Влияние изменения НВ почвы на динамику влагозапаса: 1 - влагозапас при НВ, увеличенном на 10%; 2 - влагозапас при базовом значении НВ; 3 - влагозапас при НВ, уменьшенном на 10%
На рисунках 2 и 3 приведены траектории модели при вариации наименьшей влагоемкости почвы (НВ). Данные на графиках относятся к модели яровой пшеницы, выращенной в 2005 г. на ОПХ им. Докучаева.
Траектории соответствуют базовому значению НВ, значению, измененному на ±10% для расчета динамики почвенной влаги и увеличенному на 10% для моделирования динамики биомассы органов растения.
к
с*
к
и:
Дни от момента сева
Рис. 3. Влияние величины НВ почвы на биомассу: 1 - суммарная биомасса; 2 - биомасса стеблей; 3 - биомасса зерна; 4 - биомасса листьев; 5 - биомасса корней; * - биомасса, полученная при увеличении НВ на 10%
Из рисунка 2 видно, что увеличение НВ приводит к росту влагозапасов в метровом слое почвы в течение всего вегетационного периода. Уменьшение этого параметра приводит к уменьшению влагозапасов в начальный период вегетации, которое постепенно сглаживается к концу этого периода и даже превосходит базовое значение. При увеличении НВ (рис. 3) заметно уменьшаются биомассы стеблей и зерна, биомассы корней и листьев заметно не изменяются.
Приведем результаты исследования чувствительности модели к вариации гидрофизических параметров. Исходные величины (МГ, ВЗ, НВ и ПВ) в компьютерном эксперименте варьировались как в сторону увеличения, так и уменьшения: на 20% — для МГ, ВЗ, ПВ и на 10% — для НВ. Во всех вариантах эксперимента фиксировалась величина урожая. Результаты расчета относительного коэффициента чувствительности R (16) приведены в таблице.
Первые три строки таблицы содержат средние значения урожая и соответствующие им осредненные значения коэффициента чувствительности. Последующие строки содержат среднее (ср. знач.), максимальное (макс.) и минимальное (мин.) значения R по различным вариантам. Как видно из таблицы, модель наименее чувствительна к вариации максимальной гигроскопичности. Коэффициенты чувствительности для влажности завядания и полной влагоемкости на порядок выше, чем эти же значения для МГ. Коэффициент чувстви-
тельности принимает наибольшее значение при вариации наименьшей влагоемкости (по абсолютной величине). Это означает, что к определению НВ должны быть предъявлены более высокие требования, чем к определению остальных параметров.
Значения коэффициента чувствительности модели яровой пшеницы при различных базовых уровнях урожая
Значения урожая МГ ВЗ НВ ПВ
35.18 -0.0405 -0.291 -1.19 -0.0115
28.57 0.0029 0.031 -0.20 -0.0032
23.42 0.0000 -0.017 -0.89 -0.0029
Ср. знач. -0.0125 -0.095 -0.825 -0.0059
Макс. 0.00732 0.10977 0.08337 0.07716
Мин. -0.0958 -0.5151 -5.9102 -0.4111
Значения урожая МГ ВЗ НВ ПВ
35.18 -0.0405 -0.291 -1.19 -0.0115
28.57 0.0029 0.031 -0.20 -0.0032
23.42 0.0000 -0.017 -0.89 -0.0029
Ср. знач. -0.0125 -0.095 -0.825 -0.0059
Макс. 0.00732 0.10977 0.08337 0.07716
Мин. -0.0958 -0.5151 -5.9102 -0.4111
Библиографический список
1. Пененко В. В. Методы численного моделирования атмосферных процессов. — Л., 1981.
2. Марчук Г. И. Численное решение задач динамики атмосфер и океана. — Л., 1974.
3. Полуэктов Р. А., Смоляр Э. И., Терлеев В. В., То-паж А. Г. Модели продукционного процесса сельскохозяйственных культур. — СПб., 2006.
4. Т ерлеев В. В., Полуэктов Р. А., Бакаленко Б. И. Структура информационного обеспечения модели продукционного процесса сельскохозяйственных культур // Агрофизика. — 2012 — № 2.
