On impulse control problems arising in automatization of pests control in greenhouses

Full Text

Введение

Современные инновационные решения в сельском хозяйстве направлены на то, чтобы потребители получали качественную агропродукцию, выращенную с минимальным применением химических средств. Одним из подходов к решению этой задачи является использование концепции комплексной борьбы с вредителями (Integrated Pest Management). Данная концепция подразумевает защиту агрокультур, основанную на применении комплекса мероприятий, учитывающих циклы совместного развития представленных в теплице растений, насекомых и микроорганизмов, а также их взаимодействия с питательными веществами [1, 2]. Современные подходы, направленные на снижение зависимости агрокультур от пестицидов, требуют от тепличных экосистем способности самостоятельно восстанавливаться [3]. Для создания подобных устойчивых экосистем в институте X-BIO Тюменского государственного университета ведется разработка программно-аппаратного комплекса для выработки и реализации мер профилактики и борьбы с наиболее распространенными вредителями сельскохозяйственно значимых растений, выращиваемых в закрытом грунте. Одним из модулей данного комплекса должен стать программный комплекс автоматизации экологического контроля вредителей тепличных агрокультур, использующий биопрепараты на основе энтомофагов. В основу математического обеспечения программного комплекса будут положены дифференциальные системы типа «хищник-жертва», содержащие управляющие параметры. В настоящей работе рассматривается подобная дифференциальная модель, содержащая импульсное управление, формализующее серию управляющих внесений биопрепарата-энтомофага в тепличную систему.

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

1.   Математическая модель

Рассмотрим систему

$\begin{array}{l}\frac{d}{dt}x\left(t\right)=a(1-\frac{x\left(t\right)}{K})x\left(t\right)-f\left(x\right(t\left)\right)y\left(t\right),\\ 4pt\frac{d}{dt}y\left(t\right)=cf\left(x\right(t),y(t\left)\right)y\left(t\right)-my\left(t\right),\end{array}$   (1)

где $t\in T$ ($T>0$), $x\left(t\right)$ и y(t) — значения плотности популяции вредителя и энтомофага, соответственно (в таблице 1 приведена интерпретация параметров модели). Функция  — трофическая функция Холлинга второго типа [4], имеющая вид $f\left(x\right)=\frac{bx}{d+x}.$

Восстановим значения параметров модели (1.1), опираясь на данные полевых экспериментов, проведенных сотрудниками Института исследования растений в Департаменте сельского хозяйства округа Виктория, Австралия (экспериментальные данные, а также подробное описание постановки экспериментов доступны в работе [5]), где роли тепличной агрокультуры, вредителя и энтомофага играли обыкновенный огурец Cucumis sativus, паутинный клещ Tetranychus urticae и клещ-акарифаг Tetranychus persimilis, соответственно.

 

Таблица 1

Описание параметров модели (1.1)

Параметр	Описание
$a>0$	Скорость прироста популяции вредителя
$K>0$	Плотность популяции вредителя, приводящая к полному поражению им агрокультуры
$b>0$	Максимальная скорость истребления вредителя энтомофагом
$0<c<1$	Эффективность конверсии питания энтомофага вредителем в прирост популяции энтомофага
$0<b/d<1$	Эффективность поиска вредителем энтомофага
$m>0$	Скорость естественной гибели энтомофага

 

Запишем систему (1.1) в сокращенном виде 

$\begin{array}{l}\overline{X}\left(t\right)=G\left(\overline{X}\right(t),\overline{P}),\end{array}$  (2)

где $\overline{X}\left(t\right)={\left(x\right(t),y(t\left)\right)}^T={\left(x\right(t,\overline{P}),y(t,\overline{P}\left)\right)}^T\text{},$ $G:{ℝ}^2\times M\to {ℝ}^2$ — функция, соответствующая правой части (1.1), $M\subset {ℝ}^5$ — множество допустимых значений параметров. Для решения задачи идентификации параметров системы (1.1) используем метод максимального правдоподобия. Введем функционал $\mathcal{F}\mathrm{,}$ определяющий отклонение решения модели от экспериментальных данных:

$\mathcal{F}\left(\overline{X}\right)=\sum _{t_i\in \mathrm{𝕀}}\left|\overline{X}\right(t_i,\overline{P})-{\overline{X}}_i|.$

Здесь $t_i$ — i-ая неделя эксперимента, ${\tilde{X}}_i$ — вектор экспериментальных данных по численностям популяций на неделю $t_i\mathrm{,}$ $\overline{P}$ — вектор параметров модели (см. таблицу 1), $\overline{X}(t_i,\overline{P})$ — значение  в момент времени $t_i$ при значениях параметров $\overline{P}\mathrm{,}$ полученное в результате решения системы (1.1), символом $|\cdot |$ обозначено евклидово расстояние в ${ℝ}^2\mathrm{.}$ Для идентификации параметров модели необходимо минимизировать функционал $\mathcal{F}$ относительно переменной $\overline{P}\mathrm{,}$ предварительно оценив области определения ее компонент (границы множества $M$).

Исходя из биологического смысла параметра $K\mathrm{,}$ примем его значение $K=1$ (100%). Отметим, что действительная плотность популяции вредителя может достигать значений, превышающих $\mathrm{1,}$ однако в настоящем исследовании нас интересует режим работы системы (1.1), при котором $x\left(t\right)<K$ при всех $t\in \left[\mathrm{0,}T\right],$ что соответствует сохранению тепличной культуры на рассматриваемом интервале времени (вышеупомянутые экспериментальные данные также согласуются с условием $x\left(t\right)<K$).

Оценим скорость a прироста популяции вредителя исходя из данных, собранных на экспериментальных площадках до момента введения в эксперимент клеща-акарифага(см. таблицу 2). Проинтегрируем первое уравнение системы уравнение (1.1) в отсутствие энтомофага (т. е. при $y\equiv 0$). Получим

$a=\ln \left(\frac{x\left(0\right)(K-x(1\left)\right)}{x\left(1\right)(K-x(0\left)\right)}\right).$

Будем использовать в наших вычислениях среднее значение $a=0.6159$ скорости прироста популяции вредителя.

 

Таблица 2

Экспериментальные данные по приросту численности паутинного клеща в отсутствие клеща-акарифага

Номер эксперимента	$x\left(0\right)$	$x\left(1\right)$	a
I	0.65	0.76	0.5336
II	0.36	0.56	0.8165
III	0.94	0.96	0.4265
IV	0.64	0.70	0.6897

 

Для биологически релевантной оценки границ скорости b максимального потребления добычи энтомофагом, будем опираться на данные по потреблению клещом P. persimilis яиц (см. [6–9]), личинок (см. [8, 10]) и взрослых особей (см. [8]) паутинного клеща, а также на сведения о типичной пропорции представленности вышеуказанных стадий развития паутинного клеща в популяции (см. [8]). Получим оценку $b\in \left[\mathrm{3.49,4.92}\right].$

Согласно [6], эффективность конверсии питания P. persimilis в прирост популяции имеет приблизительное значение $c\approx \mathrm{0.7.}$ В наших расчетах выберем диапазон $c\in \left[\mathrm{0.6,0.8}\right].$

Оценку значений параметра d в терминах биологически интерпретируемого параметра $b/d,$ характеризующего в рассматриваемой тепличной экосистеме эффективность поиска добычи хищником, произвести довольно затруднительно ввиду существенной зависимости этого параметра от ряда условий эксперимента (например, от температуры в тепличной системе), не представленных в формализации (1.1). В связи с этим, при восстановлении параметра d будем использовать стандартную для P. persimilis оценку $b/d\in \left[\mathrm{0.2,0.65}\right]$ (см. [5]), откуда получим $d\in \left(1.53b\mathrm{,5}b\right).$

Коэффициент смертности клеща-акарифага m полагается нами варьирующимся от 1.4 до 2.8 на основе сведений, содержащихся в работах [11, 12].

Для решения задачи поиска компонент вектора $\overline{P}\mathrm{,}$ соответствующих минимуму функционала $\mathcal{F}\mathrm{,}$ используем метод Нелдера–Мида [13]. Интегральные кривые системы (1.1) при найденном минимизирующем функционал значении ${\overline{P}}_0$ (см. таблицу 3) изображены на рисунке 1. Система (1.1) обладает единственным (неустойчивым) положением равновесия $E\mathrm{,}$ имеющим координаты $\left(\mathrm{32.67,4.99}\right)$ на фазовой плоскости и порождающим предельный цикл.

 

Таблица 3

Восстановленные по экспериментальным данным значения параметров системы (1.1)

a	b	c	d	m
0.6159	4.31	0.68	20.52	1.8

 

Рис. 1. Динамика компонент системы (1.1) при идентифицированных значениях параметров (A). Сравнение восстановленной динамики компонент x(t) (B) и y(t) (C) системы (1.1) с соответствующими экспериментальными данными (отмечены символами "*")

 

2. Задача импульсного управления

Продемонстрируем применимость модели (1.1) для автоматизации экологического контроля вредителей тепличных агрокультур. Для всякого начального условия $\left(x\right(0),y(0\left)\right)=(x_0,y_0)\in {ℝ}^2,$ уравнение (1.2) можно представить в виде

$\begin{array}{l}\overline{X}\left(t\right)={\overline{X}}_0+{\int }_a^t{G}_0\left(\overline{X}\right(s\left)\right)ds,t\in \left[\mathrm{0,}T\right],\end{array}$     (3)

 где $\overline{X}\left(0\right)={(x_0,y_0)}^T,$ $G_0\left(\overline{X}\right)=G(\overline{X},{\overline{P}}_0).$

Для системы (2.1) возможность контроля численности вредителей, основанная на внесении в тепличную систему биопрепаратов на основе энтомофагов, может быть формализована следующим образом:

$\begin{array}{l}\overline{X}\left(t\right)={\overline{X}}_0+{\int }_a^t{G}_0\left(\overline{X}\right(s\left)\right)ds+\overline{U}\left(t\right),\end{array}$   (4)

где $\overline{U}$ — импульсное управление вида ${\left(\mathrm{0,}{\sum }_{\forall k}{\chi }_{[t_k,T]}\right(t\left)y_k\right)}^T\text{},$ ${\chi }_A$ — характеристическая функция множества $A\subset \left[\mathrm{0,}T\right],$ $\left\{t_k\right\}\subset \left[\mathrm{0,}T\right]$ и $y_k>0$ — моменты времени и величины импульсных воздействий, соответственно, причем $0⩽k⩽\widehat{K}\in \mathbb{N}.$ Таким образом, различные управления  могут состоять из разного количества импульсных воздействий, общее количество которых для каждого управления не превосходит $\widehat{K}\mathrm{.}$

Определим множество y функций вида $\overline{Y}=\overline{X}+\overline{U},$ где $\overline{X}$ — функция из класса непрерывных двумерных вектор-функций. Множество y образует полное метрическое пространство относительно метрики

${\rho }_{\mathcal{Y}}({\overline{Y}}_1,{\overline{Y}}_2)=\parallel {\overline{X}}_1-{\overline{X}}_2{\parallel }_C+{\int }_0^T\left|\sum _{k:t_k^1\in [a,b]}{\chi }_{[t_k^1,b]}\right(t)y_k^1-\sum _{k:t_k^2\in [a,b]}{\chi }_{[t_k^2,b]}(t\left)y_k^2\right|dt,$  (5)

где символом $\parallel \cdot {\parallel }_C$ обозначена норма в соответствующем пространстве непрерывных функций (см. [14]). Подобная метрика была использована работе [15] при исследовании отклика нейронных систем на зрительные стимулы.

Решением задачи управления (2.2) будем считать функцию $\overline{Y}\in \mathcal{Y},$ удовлетворяющую уравнению (2.2).

Утверждение 2.1. При любом управлении $\overline{U}$ задача (2.2) имеет единственное решение. Кроме того, если некоторая последовательность управлений ${\overline{U}}^i$ сходится к управлению ${\overline{U}}^0$ в смысле сходимости ${\int }_0^T\left|{\sum }_{k:t_k^i\in \left[\mathrm{0,}T\right]}{\chi }_{[t_k^i,T]}\right(t)V_k^{i}dt-\text{}{\sum }_{k:t_k^0\in \left[\mathrm{0,}T\right]}{\chi }_{[t_k^0,T]}(t\left)V_k^0\right|dt\to 0$ при $i\to \infty ,$ то последовательность решений ${\overline{Y}}^i$ системы (2.2), соответствующих управлениям ${\overline{U}}^i\text{}\mathrm{,}$ сходится к решению ${\overline{Y}}^0\text{}\mathrm{,}$ соответствующему управлению ${\overline{U}}^0\text{}\mathrm{,}$ в метрике (2.3).

Справедливость утверждения 2.1 следует из основного результата работы [14].

Наибольший практический интерес представляет использование результатов данного исследования для решения задач экономической оптимизации экологического контроля вредителей, значимость которых неоднократно подчеркивалась в профильной литературе по направлению комплексной борьбы с вредителями (см., например, работы [16, 17], а также обзор [18] и библиографию). Утверждение 2.1 обосновывает корректность задачи импульсного управления для рассматриваемой системы «хищник-жертва», открывая возможности отыскания наиболее эффективных с экономической точки зрения управляющих стратегий, в том числе, с использованием численных методов. Подобные задачи оптимизации могут быть формализованы, например, следующим образом:

$\begin{array}{l}\overline{X}\left(t\right)={\overline{X}}_0+{\int }_a^{t}G\left(\overline{X}\right(s\left)\right)ds+\sum _{\forall k}{\chi }_{[t_k,T]}\left(t\right)y_k,t\in \left[\mathrm{0,}T\right],\\ 4pt\parallel (\overline{X},{\left(\mathrm{1,0}\right)}^T){\parallel }_C<\Theta ,\\ 4pt\sum _{\forall k}{y}_k\to \min ,\end{array}$

где $\Theta >0$ — критический уровень поражения тепличной агрокультуры вредителем, символом $(\cdot ,\cdot )$ обозначено скалярное произведение векторов.

 

Рис. 2. Шаги построения допустимого управления в рамках системы (2.2). Синие стрелки показывают направление движения по кривым фазового состояния. Черной окружностью изображено положение равновесия E системы (1.1), красной пунктирной линией изображено пороговое значение $\Theta$ поражения агрокультуры, черные пунктирные векторы отображают импульсные воздействия на экосистему теплицы. Рисунки (A), (B) и (C) отвечают первому, второму и третьему импульсному воздействию, соответственно

 

Настоящая работа не преследует цели отыскания оптимальных управляющих стратегий, однако мы приведем эвристическую схему построения стратегии импульсного управления системой «хищник-жертва», исследованной в предыдущем параграфе. Пусть $\Theta =\mathrm{0.5,}$ ${\overline{X}}_0=\left(\mathrm{0.01,0}\right).$ Фазовый портрет данной системы изображен на рисунке 2 (A). С учетом имеющегося знания о координатах положения равновесия $E=\left(\mathrm{32.67,4.99}\right)$ рассматриваемой системы, применим операцию внесения энтомофагов в тепличную систему в количестве $y_1=4.99$ в такой момент времени $t_1\mathrm{,}$что $x\left(t_1\right)\approx \mathrm{32.67.}$ В силу того, что численность популяций в тепличной экосистеме может быть оценена при помощи системы (1.1) лишь приближенно, будем полагать, что момент реализации импульсного воздействия определяется квадратом на конце черного пунктирного вектора (при этом длина вектора равна величине $y_1=\mathrm{4.99,}$ см. рисунок 2 (A)). Фазовое состояние системы, таким образом (с учетом погрешностей), переместилось в окрестность положения равновесия E (положение равновесия изображено черной окружностью, см. рисунок 2 (B)). Следующим моментом применения управления предлагается наибольшее значение $t_2\mathrm{,}$ $x\left(t_2\right)\approx \mathrm{32.67,}$ предшествующее возникновению общих точек кривой фазового состояния и линии порогового значения поражения агрокультуры (отображена красным пунктиром). С учетом неизбежных погрешностей, будем полагать, что момент импульсного воздействия определяется квадратом на конце черного пунктирного вектора (при этом длина вектора равна величине $y_2=4.99-y\left(t_2\right),$ см. рисунок 2 (B)). Фазовое состояние системы переместилось в окрестность положения равновесия  (изображено черной окружностью, см. рисунок 2 (С)). Следующим моментом применения управления предлагается наибольшее значение $t_3,$ $x\left(t_3\right)\approx 32.67$ (см. рисунок 2 (С)), предшествующее возникновению общих точек кривой фазового состояния и линии порогового значения поражения агрокультуры. Модуль  применяемого импульсного воздействия при этом равен $4.99-y\left(t_2\right).$ И так далее ... Итоговая стоимость поддержания популяции вредителя в допороговом состоянии будет пропорциональна сумме длин черных пунктирных векторов. Естественным образом, ключевую роль в снижении стоимости реализации подобных управляющих стратегий играет высокая точность оценки численности популяций вредителя и энтомофага, позволяющая, путем внесения необходимых количеств особей энтомофага в тепличный комплекс, переводить состояние системы как можно ближе к ее положению равновесия.

About the authors

Evgenii O. Burlakov

V. A. Trapeznikov Institute of Control Sciences, Russian Academy of Sciences; Derzhavin Tambov State University; Tyumen State University

Author for correspondence.
Email: eb_@bk.ru
ORCID iD: 0000-0002-7286-9456

PhD, Researcher at the V. A. Trapeznikov Institute of Control Sciences of RAS; Researcher at the Scientific and Educational Center “Fundamental Mathematical Research”, Derzhavin Tambov State University; Researcher at the Insitute X-BIO, Tyumen State University

Russian Federation, 65 Profsoyuznaya St., Moscow 117997; 33 Internatsionalnaya St., Tambov 392000; 6 Volodarskogo St., Tyumen 625003

Vitaly A. Yarema

Tyumen State University

Email: yaremavertal@mail.ru
ORCID iD: 0009-0001-7784-6720

Post-Graduate Student, Institute X-BIO

Russian Federation, 6 Volodarskogo St., Tyumen 625003

References

Supplementary files