Представим процесс движения природно-экологического (естественного) капитала в дискретном времени следующим уравнением:
$$N_{t+1} =N_t -\Delta N_t^{S-}+\Delta N_t^{N+}-\Delta N_t^{N-}+\Delta N_t^{S+}\;\;\;(1) $$Где $N_t,N_{t+1}$ - природно-экологический потенциал соответственно в момент времени $t+1,t$; $\Delta N_t^{S-}$ - уменьшение природного потенциала за счет негативной (природа разрушительной) деятельности человека (общества); $\Delta N_t^{S+}$-увеличение природного потенциала за счет позитивной (природа восстановительной) деятельности человека (общества); $\Delta N_t^{N+}$ - увеличение природного потенциала за счет естественных механизмов воспроизводства природного капитала; $\Delta N_t^{N-}$ - уменьшение природного потенциала за счет естественных причин (глобальное изменение климата, почвы и.т. д.).
Здесь $\Delta N_t^{S-}=n_t^SN_t$ и $\Delta N_t^{S+}= \Lambda(\gamma^Nknow(t)I^N (t))$ - функция трансформации инвестиций в прирост природного капитала, $\gamma^Nknow$ - характеризует уровень инновационности использования инвестиций для формирования природного капитала. При $0< \gamma^Nknow$$ < 1.0$ инвестиции используются ниже общественно необходимого уровня, а при $\gamma^Nknow$$ > 1.0$ выше. Данный коэффициент в концентрированном виде выражает качество искусственного развития природного капитала агросистемы. Для природного капитала очень важны естественные процессы роста (самовосстановления) и деградации. То есть, имеем, $\Delta N_t^{N+}=n_t^{N+} N_t$, где .$0 < \mathop n\nolimits_t^{N + } = \mathop n\nolimits_t^{regN + } + \mathop \varepsilon \nolimits_t^{N + } \le \mathop h\nolimits^{N + }$. Здесь $\mathop n\nolimits_t^{N + }$ коэффициент, характеризующий естественное самовозрастание природного капитала, $\mathop n\nolimits_t^{regN + }\;$, $\mathop \varepsilon \nolimits_t^{N + }$ регулярная и стохастическая составляющие данного коэффициента, $\mathop h\nolimits^{N + }$ -
Для $\Delta \mathop N\nolimits_t^{N-} = \mathop n\nolimits_t^{N - } \cdot \mathop N\nolimits_t$ , где $0 < \mathop n\nolimits_t^{N - } = \mathop n\nolimits_t^{regN - } + \mathop \varepsilon \nolimits_t^{N - } \le 1$. Здесь $n_t^{N-}$ коэффициент, характеризующий естественное убывание природного капитала, $\mathop n\nolimits_t^{regN - } ,\mathop \varepsilon \nolimits_t^{N - }$ регулярная и стохастическая составляющие данного коэффициента. В пределе считаем, что природный капитал может быть уничтожен за счет естественных причин, например природных катаклизмов. Следовательно, имеем:
$$k_{t+1}^N=\frac{{\mathop N\nolimits_{t + 1} }}{{\mathop N\nolimits_t }} = 1 + \frac{{\Delta \mathop N\nolimits_t^{S + } }}{{\mathop N\nolimits_t }} + \frac{{\Delta \mathop N\nolimits_t^{N + } }}{{\mathop N\nolimits_t }} - \frac{{\Delta \mathop N\nolimits_t^{S - } }}{{\mathop N\nolimits_t }} - \frac{{\Delta \mathop N\nolimits_t^{N - } }}{{\mathop N\nolimits_t }}=$$$$ = 1 + \frac{{\Lambda (\mathop \gamma \nolimits_t^{Nknow} \mathop I\nolimits_t^N )}}{{\mathop N\nolimits_t }} + \mathop n\nolimits_t^{N + } - \mathop n\nolimits_t^{S - } - \mathop n\nolimits_t^{N - }\;\;\;(2)$$Из выражения (2) видно, что развитие природно-экологического капитала определяется естественными природными процессами, формирующими естественный прирост и убыль природного капитала $(n_t^{N+}, n_t^{N-})$, социально политическими и экономическими процессами, детерминирующими негативное антропогенное воздействие на природный капитал $(n_t^{S-})$, уровнем инвестиций в природоохранную деятельность $(I_t^N)$ и уровнем инновационности использования инвестиций в развитие природного капитала $(\gamma_t^{Nknow})$. Основное дискретное динамическое соотношение:
$$\mathop N\nolimits_{t + 1} = \mathop k\nolimits_{t + 1}^N \cdot \mathop N\nolimits_t\;\;\;(3)$$Перейдем от дискретного к не прерывному времени. Вычитая из обеих частей уравнения (3) $N_t$ получим следующее уравнение:
$$\mathop N\nolimits_{t + 1} - \mathop N\nolimits_t = \mathop {\Delta \mathop N\nolimits_t = (k}\nolimits_{t + 1}^N - 1) \cdot \mathop N\nolimits_t$$$$\Delta \mathop N\nolimits_t = \mathop {\tilde k}\nolimits_{t + 1}^N \cdot \mathop N\nolimits_t\;\;\;(4)$$Переходя к бесконечно малым приращениям, получаем следующее дифференциальное уравнение:
$$\frac{{dN}}{{dt}} = \mathop {\tilde k}\nolimits_{}^N (t) \cdot \mathop N\nolimits^{} (t) = \Lambda (\mathop \gamma \nolimits^{Nknow} (t),\mathop I\nolimits^N (t)) + \mathop N\nolimits^{} (t)(\mathop n\nolimits^{N + } (t) - \mathop n\nolimits^{S - } (t) - \mathop n\nolimits^{N - } (t))\;\;\;\;(5)$$Получим удельное динамическое уравнение для человеческого капита-ла в дискретном времени:
$\mathop n\nolimits_{t + 1} = (\frac{{\mathop k\nolimits_{t + 1}^N }}{{(1 + \mathop \omega \nolimits_t )}} - 1) \cdot \mathop n\nolimits_t$ , где $n_{t+1}$ и $n_t$ природный капитал на душу населения соответственно в $t+1$ и $t$ периоды времени.
$\Delta \mathop n\nolimits_t = (\frac{{\mathop k\nolimits_{t + 1}^N }}{{(1 + \mathop \omega \nolimits_t )}} - 1) \cdot \mathop n\nolimits_t$ и переходя к малым приращениям, получаем удельное дифференциальное уравнение:
Рассмотрим соотношение (5), сформулировав некоторые экономически и экологически правдоподобные гипотезы, задающие соответствующие тренды членов правой части данного дифференциального уравнения. Обратимся сначала к функции, являющейся переменным коэффициентом при уровне природного капитала АСЭЭС $N(t)$ в момент $t$:
$$W(t) = n^{N + } (t) - n^{S - } (t) - n^{N - } (t)\;\;\;(7)$$Здесь представлены три составляющих $n^{N+}(t)$ и $n^{N-}(t)$ характеризуют естественные процессы повышения (снижения) размера и качества природного капитала АСЭЭС, а $n_t^{S-}(t)$ описывает его разрушение в результате человеческой деятельности. Прежде всего, необходимо определить само понятие «природный капитал аграрного производства». Приведем здесь только предварительные соображения по данному вопросу, которые достаточны для демонстрации возможностей сформулированных выше моделей динамики природного капитала агросферы. И так, в качестве природного капитала нами рассматривается подсистема природной системы, которая включается в процесс аграрного производства в качестве предмета или орудия труда, то есть это те природные объекты и силы, без которых невозможно в современных социально-экономических условиях производство сельскохозяйственной продукции. Естественной основой аграрного производства является территория, обладающая определенными топологическими и климатическими свойства (форма, рельеф, площадь, типы почв, наличие водоемов, годовое количество солнечной энергии, температурный режим, уровень влажности и т.п). Неравномерность распределения этих свойств по территории приводит к формированию ареалов, где имеется такой набор этих свойств, что естественные природные силы возможно и целесообразно использовать в данных кокретно-исторических условиях для сельскохозяйственного производства. Именно эти ареалы территории - есть естественная основа природного капитала АСЭЭС, в качестве которого выступает, прежде всего, земля сельскохозяйственного назначения. Основными характеристиками, отражающими производительные возможности природного капитала АСЭЭС, являются размер земельных угодий (площадь), плодородие земли, биоклиматический потенциал (БКП), уровень аграрной освоенности территории (отношение площади земель аграрного назначения к площади территории), пространственно-топологические характеристики (форма, расположение относительно центров социальной и экономической активности) и инфраструктурная обеспеченность. Поэтому при оценке возможных трендов $n^{N+}(t)$ и $n^{N-}(t)$ на ближайшие 20-40 для АСЭЭС России мы будем анализировать возможную динамику именно этих составляющих природного капитала. Выше мы ввели соотношение $n^{N - } (t) = n^{regN - } (t) + \varepsilon ^{N - } (t)$ , если рассматривать регулярную составляющую, то за сценарный период (20-40 лет) трудно предполагать наличие отрицательной естественной динамики. Сокращение площади сельскохозяйственных угодий России за счет естественных причин (поднятие уровня мирового океана, опустынивание и.т. д.), маловероятно, по крайней мере, общепризнанных прогнозов такого плана нет. БКП также не имеет явных тенденций к изменению (глобальное изменение климата если даже и происходит, то в рассматриваемом временном масштабе оно, по нашему мнению, ничтожно мало), не предвидится естественное отрицательное изменение уровня естественного плодородия почв и соответствующих пространственно-топологических характеристик [130,131]. Следовательно, можно принять $n^{regN - } (t) \approx 0$, случайная составляющая связана в основном с климатическими факторами, а они, как показывают многолетние наблюдения, являются нормально распределенными случайными величинами с нулевым математическим ожиданием и не очень большой на периоде 20-40 лет дисперсией, поэтому полагаем $\varepsilon ^{N - } (t) \approx 0 \to n^{N - } (t) \approx 0$ Что касается процессов естественного роста природного капитала $n^{N+}(t)$, то проблема здесь достаточно сложна и общепризнанных подходов системного интегрального характера, по нашему мнению нет, мы все же постараемся их учесть в рамках предлагаемых моделей. Но вполне очевидно, что в модель можно включать, если это необходимо, эффекты $n^{N-}(t)$ и $n^{N+}(t)$ как в дискретном, так и непрерывном вариантах.
Рассмотрим, как формируется отрицательный виртуальный прирост природного капитала АСЭЭС за счет уменьшения продуктивности земли. Пусть в момент $t$ $y_t$ – продуктивность единицы площади земли и $S_t$ площадь сельскохозяйственных угодий (полагаем $S_t=const)$, тогда $\mathop P\nolimits_{t} = \mathop y\nolimits_t \mathop S\nolimits_t$ , где $P_t$ объем сельскохозяйственной продукции, аналогично имеем $\mathop P\nolimits_{t + 1} = \mathop y\nolimits_{t + 1} \mathop { \cdot S}\nolimits_t$ При этом имеет место следующее соотношение $\mathop y\nolimits_{t + 1} = \mathop y\nolimits_t \cdot \mathop a\nolimits_t\;\;\;$ $\mathop a\nolimits^{pr} < \mathop a\nolimits_t \le 1$ –коэффициент характеризующий темп падения продуктивности земли момент времени $t$, тогда $\mathop P\nolimits_{t + 1} = \mathop y\nolimits_t \mathop { \cdot \mathop a\nolimits_t \cdot S}\nolimits_t$ отсюда:
$\mathop {\Delta P}\nolimits_t = \mathop S\nolimits_t \cdot \mathop y\nolimits_{t - 1} \cdot (1 - \mathop a\nolimits_t )$,следовательно, имеем $\mathop S\nolimits_{t + 1}^V \cdot \mathop y\nolimits_t \cdot \mathop a\nolimits_t = \mathop P\nolimits_{t + 1} + \mathop {\Delta P}\nolimits_t$ отсюда $\mathop S\nolimits_{t + 1}^V = \frac{{\mathop S\nolimits_t }}{{\mathop a\nolimits_t }}\;\;\;$, где $S_{t+1}^v$ где виртуальная площадь сельхозугодий, которую нужно было бы иметь, чтобы получить объем сельхозпродукции равный $P_t$.Тогда виртуальный прирост природного капитала:
$$\Delta \mathop N\nolimits_t^{S - } = \mathop {\Delta S}\nolimits_t = \mathop S\nolimits_t - \frac{{\mathop S\nolimits_t }}{{\mathop a\nolimits_t }} = \mathop S\nolimits_t \cdot (1 - \frac{1}{{\mathop a\nolimits_t }}) = \mathop N\nolimits_t (1 - \frac{1}{{\mathop a\nolimits_t }}) \le 0\;\;\;\;(8)$$Величина $a_t$, определяющая продуктивность земли в любой момент времени, зависит от экологической нагрузки, которая оказывает экономика на природный капитал АСЭЭС. В качестве измерителя нагрузки возьмем относительную экологическую нагрузку как отношение текущей нагрузки на максимально допустимую, то есть такую нагрузку, при которой уровень продуктивности земли таков, что начинается резкое снижение качества жизни в связи с невозможностью обеспечить население продуктами питания в необходимом количестве и качестве. Определение придельной экологической нагрузки рассмотрим ниже. Ясно, что относительная экологическая нагрузка должна представлять собой монотонно возрастающую функции времени с интервалом значений $[0.0,1.0]$то есть, имеем: $\mathop D\nolimits^{otn} (t) \in [0.0,1.0],{\rm t} \in {\rm [0}{\rm , + }\infty {\rm )}$.
Коэффициент темпа падения продуктивности земли $a_t$ из естественных соображений представляет собой монотонно убывающую функцию $D^{otn}$ на интервале $(a^{pr},1.0]$. При этом в отсутствии нагрузки он должен принимать значение равное 1 , при придельной нагрузке значение $a^{pr}$, определяющее максимальный темп падения продуктивности земли, при котором фактически начинается структурно-функциональное разрушение национальной АСЭЭС. Предельное значение $a^{pr}=0.5$ (отрицательный прирост природного капитала при этом равен его первоначальному значению). Можно предложить следующие классы моделей падения продуктивности земли от уровня относительной экологической загрузки:
$\mathop a\nolimits_t = - \mathop a\nolimits^{pr} \mathop D\nolimits_{otn}^{(m)} + 1\;\;\;$ где параметр $m$ задает скорость и характер падения $a_t$. При $m=1$ имеем линейное падение с постоянной скоростью ( а), $при m=1.5$ нелинейное падение с увеличивающейся скоростью (б) и при $\;\;\;m=0.5$ нелинейное падение с уменьшающейся скоростью (в).
import warnings
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import rand
warnings.simplefilter('ignore')
Dotn = np.arange(0, 1, 0.01)
apr=[0.1,0.25,0.4,0.5]
Colat=['g','c','y','coral']
L=[r'$a^{pr}=0.1$',r'$a^{pr}=0.25$',r'$a^{pr}=0.4$',r'$a^{pr}=0.1$']
#plt.figure(figsize=(12, 8))
fig, axs = plt.subplots(1, 3, figsize=(24, 8))
m=1.0
for i in range(4):
at=-(apr[i]*Dotn**m)+1
axs[0].plot(Dotn,at,color=Colat[i], linewidth=2.0)
axs[0].set_xlabel(r'$D^{otn}$',fontsize=16)
axs[0].set_ylabel(r'$a_t$',fontsize=16)
axs[0].set_title(r'$a_t=-a^{pr}D_{otn}^{(m)}+1$', fontsize=20)
axs[0].legend(L, fontsize=14)
axs[0].grid(True)
axs[0].text(0.0,0.75,'m=1.0', fontsize=20)
m=1.5
#plt.figure(figsize=(12, 8))
for i in range(4):
at=-(apr[i]*Dotn**m)+1
axs[1].plot(Dotn,at,color=Colat[i], linewidth=2.0)
axs[1].set_xlabel(r'$D^{otn}$',fontsize=16)
axs[1].set_ylabel(r'$a_t$',fontsize=16)
axs[1].set_title(r'$a_t=-a^{pr}D_{otn}^{(m)}+1$', fontsize=20)
axs[1].legend(L, fontsize=14)
axs[1].grid(True)
axs[1].text(0.0,0.75,'m=1.5', fontsize=20)
m=0.5
#plt.figure(figsize=(12, 8))
for i in range(4):
at=-(apr[i]*Dotn**m)+1
axs[2].plot(Dotn,at,color=Colat[i], linewidth=2.0)
axs[2].set_xlabel(r'$D^{otn}$',fontsize=16)
axs[2].set_ylabel(r'$a_t$',fontsize=16)
axs[2].set_title(r'$a_t=-a^{pr}D_{otn}^{(m)}+1$', fontsize=20)
axs[2].legend(L, fontsize=14)
axs[2].grid(True)
axs[2].text(0.0,0.75,'m=0.5', fontsize=20)
fig.tight_layout()
plt.show()
Перейдем к оценке относительной экологической нагрузки. В качестве измерителя предлагаем использовать следующее выражение: $\frac{{D(P(t))}}{{\mathop D\nolimits_{ПН} }}$ , где P(t) - функции описывающая динамику ВВП ; $D(P(t))$ – антропогенная нагрузка, разрушающая природный потенциал, как функция ВВП; $D_{ПН}$ – предельно допустимая нагрузка, при которой природный потенциал АСЭЭС разрушается полностью. Для оценки этой функции для России на период до 50 лет воспользуемся результатами глобального моделирования мирового развития [Римский клуб]. Для иллюстрации предлагаемых подходов выбираем два сценария пессимистический (сценарий 7) (рис.) и оптимистический (сценарий 9) (рис. ):


Сценарий 9 предполагает в перспективе постепенную стабилизацию населения, снижение темпов роста экономики, повышение ее инновационности и экологичности. В сценарии 7 предполагается сохранение, и даже усиление существующих негативных тенденций развития. На графике сценария 7 приблизительно на 2053 год приходится максимум экологической нагрузки, которая равна 106 относительных единиц. Будем считать, что это предельная экологическая нагрузка человечества: $D_{ПН}^{WORLD}=106 \;o.e.$
Для определения значения предельной экологической нагрузки для России определим прогноз доли ВВП России в мировом ВВП на 2053 год. По данным ЦРУ США ВВП России по паритету покупательной способности в 2010 году составил 2.99 процента от мирового. Если считать, что мировой ВВП будет расти в среднем 2 процента в год, а ВВП России в среднем 4 процента в год, то доля России в мировом ВВП к 2053 году может составить порядка 6.9 процента. Следовательно, можно принять $D_{ПН}=0.069D_{ПН}^{WORLD}=7.3\; o.e.$
Для построения функции $D(P(t))$ сначала сконструируем функцию $P(t)$ прогноз изменения ВВП России на период 2010-2053 гг. Будем исходить из оптимистического предположения, что среднегодовой темп роста составит 4 процента, тогда $P_0=2229$ млн. долларов США (фактически 2010 г.), а в конце сценарного срока 2053 год $P_{43}=12038$ млн. долларов США. Наиболее вероятной динамической моделью ВВП России на данный период лет мы считаем логистическую функцию, хорошо описывающую быстро растущие экономики в период перехода мировой экономике к повышательной фазе длинноволнового К-цикла (переход к 6 технологическому укладу). Следовательно:
$$P(t) = \frac{{12038 \cdot 2229 \cdot \mathop e\nolimits^{0.2 \cdot t} }}{{12038 + 2229 \cdot \mathop {(e}\nolimits^{0.2 \cdot t} - 1\mathop )\nolimits_{}^{} }}\;\;\;(9)$$Что касается экологической нагрузки, то, при принятом выше, сценарии динамики ВВП России целесообразно использовать откорректированный по ВВП сценарий 7 глобального развития, где начальное и конечное значения равны $D_0=1.65$ о.е. ,а $D_{43}=2.8$ о.е. :
$$\mathop D\nolimits^ + (t) = \frac{{2.8 \cdot 1.65 \cdot \mathop e\nolimits^{0.075 \cdot t} }}{{2.8 + 1.65 \cdot \mathop {(e}\nolimits^{0.075 \cdot t} - 1\mathop )\nolimits_{}^{} }}\;\;\;(10)$$import warnings
import matplotlib.pyplot as plt
import numpy as np
warnings.simplefilter('ignore')
t= np.arange(0, 43, 1)
#plt.figure(figsize=(12, 8))
fig, axs = plt.subplots(1, 2, figsize=(24, 10))
Pt=(12038*2229*np.exp(0.2*t))/(12038+2229*(np.exp(0.2*t)+1))
axs[0].plot(t,Pt,color='black', linewidth=3.0)
axs[0].set_xlabel(r'$t$',fontsize=22)
axs[0].set_ylabel(r'$P$ млн.дол.США',fontsize=22)
axs[0].set_title(r' Прогноз ВВП России 2010-2053 гг. ', fontsize=25)
axs[0].grid(True)
Dt=(2.8*1.65*np.exp(0.075*t))/(2.8+1.65*(np.exp(0.075*t)+1))
axs[1].plot(t,Dt,color='red', linewidth=2.0)
axs[1].set_xlabel(r'$t$',fontsize=22)
axs[1].set_ylabel(r'$D$ o.e.',fontsize=22)
axs[1].set_title(r'Прогноз экологической нагрузки экономики России 2010-2053 гг.', fontsize=25)
axs[1].grid(True)
fig.tight_layout()
plt.show()
Экономический смысл заключается в следующем, первые 30 лет используются относительно экологоемкие факторы роста, быстро увеличивающие экологическую нагрузку, а затем более инновационные и экологически чистые способы повышения ВВП.
Зависимость $D(P)$, построенная на основе выше вышеприведенных соотношений имеет вид:
import warnings
import matplotlib.pyplot as plt
import numpy as np
warnings.simplefilter('ignore')
t= np.arange(0, 20, 1)
Pt=(12038*2229*np.exp(0.2*t))/(12038+2229*(np.exp(0.2*t)+1))
Dt=(2.8*1.65*np.exp(0.075*t))/(2.8+1.65*(np.exp(0.075*t)+1))
mod = np.polyfit(Pt,Dt,1)
Dmod= mod[0]*Pt +mod[1]
plt.figure(figsize=(12, 8))
plt.plot(Pt,Dmod,color='blue', linewidth=2.0)
plt.plot(Pt,Dt,color='r', linewidth=2.0)
plt.xlabel(r'$P_t$ млн.дол.США',fontsize=18)
plt.ylabel(r'$D_t ,D_{mod}$ o.e.',fontsize=18)
plt.title(r'Зависимость Dt и Dmod от величины ВВП России(2010-2030 гг.)', fontsize=18)
plt.grid(True)
plt.legend([r'$D_t$',r'$D_{mod}=a_1P_t+a_0$'], fontsize=18)
plt.show()
t= np.arange(0, 43, 1)
Pt=(12038*2229*np.exp(0.2*t))/(12038+2229*(np.exp(0.2*t)+1))
Dt=(2.8*1.65*np.exp(0.075*t))/(2.8+1.65*(np.exp(0.075*t)+1))
Dmod= mod[0]*Pt +mod[1]
plt.figure(figsize=(12, 8))
plt.plot(Pt,Dt,color='r', linewidth=2.0)
plt.plot(Pt,Dmod,color='blue', linewidth=2.0)
plt.xlabel(r'$P_t$ млн.дол.США',fontsize=18)
plt.ylabel(r'$D_t,D_{mod}$ o.e.',fontsize=18)
plt.title(r'Зависимость Dt,Dmod от величины ВВП России (2010-2053 гг.)', fontsize=18)
plt.grid(True)
plt.legend([r'$D_t$',r'$D_{mod}$ получена на [0,10000]'], fontsize=18)
plt.show()
mod
array([9.83460059e-05, 6.12207594e-01])
Из графиков видно, что зависимость экологической нагрузки от величины ВВП есть возрастающая функция, которая имеет точку P(2030)=10000 смены тенденций, т.е от 2010 г. до 2030 г. темп роста $Dt$ практически постоянный так как D(P) очень хорошо описывается линейной зависимостью и $\frac{{dD}}{{dP}} \approx {\rm 9}{\rm .83460059e - 05}$,но в дальнейшем скорость роста экологической нагрузки резко возрастает, хотя и остается в сценарный период досточно далеко от $D_{ПН}$
Интересно сравнить экологические нагрузки, которую может оказать экономика России при реализации сценария 7 (неустойчивое разрушающее развитие) и сценария 9 (устойчивое развитие). Для этого по данным глобальной модели Медоуза, скорректированной по ВВП для России построим модель для разрушительного сценария $D^- (t)$. Эта модель строится из предположения экспоненциального роста экологической нагрузки, параметры функции оценивались на основе аппроксимации откорректированных для России данных глобальной модели Медоуза.
Поделив функции $D^+ (t)$ и $D^- (t)$, на предельное значение, можно получить соответствующие относительные функции экологической нагрузки.
Для анализа поведения природного капитала в зоне разрушения введем в рассмотрение сценарий «Апокалипсис», в котором предельный уровень достигается через 25 лет. В результате получаем три относительных функции экологической нагрузки:
$\mathop D\nolimits_{ОТН}^{ - - } (t) = 0.229 \cdot \exp (0.06 \cdot t)$ - предельно пессимистический сценарий;
$\mathop D\nolimits_{ОТН}^ - (t) = 0.2289 \cdot \exp (0.0344 \cdot t)$ - пессимистический сценарий;
$\mathop D\nolimits_{ОТН}^ + (t) = \frac{{0.6329 \cdot \exp (0.075 \cdot t)}}{{2.8 + 1.65 \cdot ( - 1 + \exp (0.075 \cdot t))}}$ - оптимистический сценапий.
Графически эти сценарии представлены ниже
import warnings
import matplotlib.pyplot as plt
import numpy as np
warnings.simplefilter('ignore')
t= np.arange(0, 43, 1)
D1t=0.229*np.exp(0.06*t)
D2t=0.2289*np.exp(0.0344*t)
D3t=0.6329*np.exp(0.075*t)/(2.8+1.65*(np.exp(0.075*t)-1))
plt.figure(figsize=(12, 8))
plt.plot(t,D1t,color='r', linewidth=2.0)
plt.plot(t,D2t,color='y', linewidth=2.0)
plt.plot(t,D3t,color='green', linewidth=2.0)
plt.xlabel(r'$t$ годы',fontsize=18)
plt.ylabel(r'$D_{ОТН}^{--}(t) ,D_{ОТН}^{-}(t),D_{ОТН}^{+}(t)$ o.e.',fontsize=18)
plt.title(r'Сценарии по экологической нагрузке в России (2010-2053 гг.)', fontsize=18)
plt.grid(True)
plt.legend([r'$D_{ОТН}^{--}(t)$ - Апокалипсис',r'$D_{ОТН}^{-}(t)$ – пессимистический сценарий',
r'$D_{ОТН}^{+}(t)$ – оптимистический сценарий'],fontsize=18)
plt.hlines(1.0, 0, 43,
color = 'r',
linewidth = 2,
linestyle = '--')
plt.vlines(25, 0, 1.0,
color = 'r',
linewidth = 1,
linestyle = '--')
plt.vlines(42, 0, 1.0,
color = 'r',
linewidth = 1,
linestyle = '--')
scatter1 = plt.scatter([25,42],[1.0,1.0],s=20* 10,c='r')
plt.annotate ("Точки выхода в зону недопустимой\n экологической нагрузки (25,1.0) (42,1.0)",
xy=(25, 1.05), xycoords='data',
xytext=(11.8, 1.7), textcoords='data',
arrowprops=dict(arrowstyle="->", connectionstyle="arc3",color='r'),
fontsize=14)
plt.annotate (" ",
xy=(42, 1.05), xycoords='data',
xytext=(15, 1.7), textcoords='data',
arrowprops=dict(arrowstyle="->", connectionstyle="arc3",color='r'))
plt.annotate ("(25,0.33)",
xy=(25, 0.33), xycoords='data',
xytext=(27, 0.15), textcoords='data',
arrowprops=dict(arrowstyle="->", connectionstyle="arc3",color='g'),fontsize=14)
plt.annotate ("(42,0.37)",
xy=(42, 0.37), xycoords='data',
xytext=(38, 0.55), textcoords='data',
arrowprops=dict(arrowstyle="->", connectionstyle="arc3",color='g'),fontsize=14)
plt.annotate ("(25,0.54)",
xy=(25, 0.547), xycoords='data',
xytext=(26, 0.83), textcoords='data',
arrowprops=dict(arrowstyle="->", connectionstyle="arc3",color='y'),fontsize=14)
scatter1 = plt.scatter([25,42],[0.334,0.37],s=8* 10,c='g')
scatter1 = plt.scatter(25,0.54,s=8* 10,c='y')
plt.show()
Из (8) имеем:
$$\mathop k\nolimits_t^{экология} = (1 - \frac{1}{{\mathop a\nolimits_t (\mathop D\nolimits_{OTH} (t))}}) \le 0\;\;\;(11)$$При этом $\mathop k\nolimits_t^{экология} = \mathop n\nolimits^{N + } (t) - \mathop n\nolimits^{S - } (t) - \mathop n\nolimits^{N - } (t)$, где $\mathop n\nolimits^{N - } (t) \approx 0,\;\;\;\mathop n\nolimits^{N + } (t) \approx 0$
Для того что бы учесть естественную восстанавливающую силу природы введем коэффициент демпфирования вредного влияния экологической нагрузки на природный капитал АСЭЭС $(n^{N+}(t)>0 )$ : $0<k_{dem}<= 1$
При $k_{dem}=1$ восстанавливающая сила природы равна 0, а при $\mathop k\nolimits_{dem} \to 0$ она стремиться к некоторому предельному уровню. Вопрос состоит в том, как включить этот параметр в модель динамики и на основе, каких предположений оценивать:
Коэффициент демпфирования должен уменьшать падение уровня продуктивности земли в любой точке интервала $t∈[0,t_{pr}]$, где $t_{pr}$ - конечная точка сценария.
При уменьшении данного коэффициента функция $k_t^{экология}$ должна смещаться вверх, но при этом оставаться отрицательной, что должно обеспечивать не возрастающий характер функции $N(t)$ (природный капитал не может самовозрастать под антропогенной экологической нагрузкой)
В пространстве параметров должна существовать область (в предельном случае точка), в которой модель динамики обеспечивает выход $N(t)$ в заданные временные диапазоны в зону разрушения природного капитала при соответствующих сценариях (сценарий 1 25-30 годы, сценарий 2 39-43 годы)
Должны существовать $a_t, k_{dem}$ обеспечивающие в любой точке интервала $t∈[0,t_{pr}]$ выполнение условия:
Включение $k_{dem}$ в $W(t)=k_t^{экология}$, при котором выполняются приведенные выше условия, может быть выполнено следующими способами:
$$a.\;\;\;W(t) = \mathop k\nolimits_t^{экология} = (1 - \frac{{\mathop k\nolimits_{dem} }}{{\mathop a\nolimits_t (\mathop D\nolimits_{OTH}^m (t))}})\;\;\;\;(13)$$$$б.\;\;\;\;W(t) = \mathop k\nolimits_t^{экология} = (1 - \frac{1}{{ - \mathop a\nolimits^{pr} \cdot \mathop k\nolimits_{dem} \cdot \mathop D\nolimits_{OTH}^m (t) + 1}})\;\;\;\;(14)$$Для случая (а) можно записать:
$$W(t) = \mathop k\nolimits_t^{экология} = \mathop n\nolimits^{N + } (t) - \mathop n\nolimits^{S - } (t)\;\;\;(15)$$Отсюда
$$\mathop n\nolimits^{N + } (t) = (\frac{1}{{\mathop a\nolimits_t (\mathop D\nolimits_{OTH}^m (t))}}) \cdot (1 + \mathop k\nolimits_{dem} )\;\;\;(16)$$то есть $0\mathop { \le n}\nolimits^{N + } (t) \le \frac{1}{{\mathop a\nolimits_t (\mathop D\nolimits_{OTH}^m (t))}}$
Или для случая (б):
$$\mathop n\nolimits^{N + } (t) = ((\frac{1}{{ - \mathop a\nolimits^{pr} \cdot \mathop D\nolimits_{OTH}^m (t)}}) + \frac{1}{{ - \mathop a\nolimits^{pr} \cdot \mathop k\nolimits_{dem} \cdot \mathop D\nolimits_{OTH}^m (t) + 1}})\;\;\;\;(17)$$Таким образом, имея значение коэффициента демпфирования, можно оценивать положительную динамику влияния процессов самоочищения и самовосстановления природного капитала (см. ниже)):
import warnings
import matplotlib.pyplot as plt
import numpy as np
warnings.simplefilter('ignore')
m=1.45
apr=0.5
kdem=[0.5,0.7,0.9]
Dotn= np.arange(0.1, 1, 0.01)
ColnN=['g','c','y']
L=[r'$k_{dem}=0.5$',r'$k_{dem}=0.7$',r'$k_{dem}=0.9$']
fig, axs = plt.subplots(2, 1, figsize=(12, 16))
for i in range(3):
at=-apr*Dotn**m+1
k_ek =(kdem[i]/at)
axs[0].plot(Dotn,k_ek,color=ColnN[i], linewidth=3.0)
axs[0].set_xlabel(r'$D_{otn}$',fontsize=20)
axs[0].set_ylabel(r'$k_t^{экология}$',fontsize=20)
axs[0].set_title(r'$k_t^{экология}$ вариант а)', fontsize=20)
axs[0].legend(L, fontsize=20)
axs[0].grid(True)
axs[0].text(0.2,1.25,'m=1.45', fontsize=20)
axs[0].text(0.2,1.15,r'$a^{pr}=0.5$', fontsize=20)
for i in range(3):
k_ek =(1/(-apr*kdem[i]*Dotn**m+1))
axs[1].plot(Dotn,k_ek,color=ColnN[i], linewidth=3.0)
axs[1].set_xlabel(r'$D_{otn}$',fontsize=20)
axs[1].set_ylabel(r'$k_t^{экология}$',fontsize=20)
axs[1].set_title(r'$k_t^{экология}$ вариант б)', fontsize=20)
axs[1].legend(L, fontsize=20)
axs[1].grid(True)
axs[1].text(0.2,1.45,'m=1.45', fontsize=20)
axs[1].text(0.2,1.35,r'$a^{pr}=0.5$', fontsize=20)
fig.tight_layout()
plt.show()
Варианты (а) и (б) демонстрируют две разные «стратегии» природы во взаимоотношении с человеком. При слабой величине восстановительной силы природы $(0.7<=k_{dem}<=1.0)$ в обоих вариантах $k^{экология}(t)$ вблизи предельной нагрузки имеет близкие значения, то есть в варианте (б) реализуется в данной области $k_{dem}$ значительно более высокий темп роста $k^{экология}(t)$ при повышении экологической нагрузки. В интервале $0.0<=k_{dem} <=0.5$ различия в темпе роста $k^{экология}(t)$ сначала снижаются, а при $k_{dem} <=0.1$ варианты (а) и (б) начинают показывать практически одинаковые темпы роста $k^{экология}(t)$.
Но при этом быстро растет различие в начальном значении $k^{экология}(t)$ в результате для $k_{dem}<=0.1$ разница в значениях $k^{экология}(t)$ около предельной нагрузки достигает различия несколько раз в пользу варианта (б). Таким образом, вариант (б) обеспечивает имитацию более мощного позитивного воздействия природы особенно в условиях высокой экологической нагрузки.
Если говорить об условиях сформулированных выше, то для варианта (б) условия 1 и 2 выполнятся вполне очевидно, для варианта (а) условие 1 выполняется в силу того, что уменьшение числителя в дроби выражения (13) тождественно увеличению знаменателя в этом же выражении, а это эквивалентно росту продуктивности земли. Условие 2 для варианта (а) требует специального анализа. Построим в системе координат $(k_{dem}, a_t)$ изокванты для которых выполняются условия: $\mathop k\nolimits_t^{экология} = (1 - \frac{{\mathop k\nolimits_{dem} }}{{\mathop a\nolimits_t (\mathop D\nolimits_{OTH} (t))}}) = const$
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
kdem=np.arange(0, 2, 0.01)
at=np.arange(0.01, 1.0, 0.01)
kdem, at = np.meshgrid(kdem,at)
kt=1-(kdem/at)
lev = [-1.0,-0.75,-0.5,-0.25,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,1.0]
plt.figure(figsize=(14, 10))
cs=plt.contour(at,kdem, kt, levels = lev, colors= ['r','r','r','r','black',
'g','g','g','g','g','g','g','g','g'],linewidths=0.75)
cs.clabel(colors='black',fmt= r'$k_t=%.2f$',fontsize=18)
plt.xlabel(r'$a_{t}$',fontsize=20)
plt.ylabel(r'$k_{dem}$',fontsize=20)
plt.title(r'Изокванты $k_t^{экология}=(1-k_{dem}/a_t)=const$', fontsize=20)
plt.grid(True)
plt.hlines(1.0, 0, 1.0,
color = 'r',
linewidth = 1,
linestyle = '--')
plt.vlines(0.5, 0, 2.0,
color = 'r',
linewidth = 1,
linestyle = '--')
plt.vlines(0.9, 0, 2.0,
color = 'r',
linewidth = 1,
linestyle = '--')
plt.hlines(0.89, 0, 1.0,
color = 'g',
linewidth = 1,
linestyle = '--')
plt.text(0.15,2.0,'РАЗРУШЕНИЕ', fontsize=20)
plt.text(0.6,2.0,'ДЕГРАДАЦИЯ', fontsize=20)
plt.text(0.93,2.0,'РОСТ', fontsize=20)
plt.annotate (r"$k_{dem}=0.9$",
xy=(0.1, 0.9), xycoords='data',
xytext=(0.01, 0.75), textcoords='data',
arrowprops=dict(arrowstyle="->", connectionstyle="arc3",color='black'),fontsize=18)
scatter1 = plt.scatter(0.98,0.89,s=11* 10,c='g')
plt.vlines(0.98, 0, 0.9,
color = 'g',
linewidth = 1,
linestyle = '--')
plt.show()
Из графика видно, что только при $k_{dem}=1.0$ во всем диапазоне изменения $a_t$ обеспечивается отрицательные значения $k_t^{экология}$ - движение по изоквантам от нулевой до больших отрицательных. При уменьшении $k_{dem}$ появляются области значений $a_t$ , где $k_t^{экология}$ принимает положительные значения. Поэтому целесообразно выбирать уровень начальной антропогенной нагрузки отличным от нуля $(a_0=A<1)$ и тогда находя точку пересечения этого значения с нулевой изоквантой определяем нижнюю границу интервала значений $k_{dem}$ в котором выполняется условие 2.
Учитывая выше приведенные рассуждения и воспользовавшись выражением (5) получаем основное уравнение динамики природного капитала (без природо-восстановительной деятельности человека):
$$\frac{{dN}}{{dt}} = \mathop W\nolimits_i (t,\mathop k\nolimits_{dem} ,\mathop a\nolimits^{pr} ,m) \cdot N(t),\;\;\;i = 1,3\;\;\; индекс\;\;\;сценария\;\;\;\;(17)$$Задавая различные значения коэффициента демпфирования, а также $a^{pr}$ и $m$ можно изучать поведение системы при разных сценариях развития экологической ситуации в АСЭЭС и разных гипотезах влияния экологии на динамику природного капитала. Используемые в динамической модели, как константы, параметры, $k_{dem},a^{pr},m$ безусловно, являются сложными функциями, изменения которых во времени связаны с фундаментальными законами природы, формулируемыми экологией, биологией, почвоведением, агрохимией и другими естественными науками. Здесь еще много неизвестного, но ряд уже полученных результатов позволяет формулировать гипотезы о динамике адаптационных свойств природы. При этом изучение влияния предложенных параметров на динамику природного капитала позволяет проверять уже имеющиеся и выдвигать дополнительные гипотезы о механизмах взаимовлиянии человека и природы. Располагая основной динамической моделью, можно перейти к моделированию динамики природного капитала. Но предварительно введем формальное определение понятий разрушение природного капитала АСЭЭС и точки разрушения.
Природный капитал разрушился в момент $t_{раз}$, если выполняется условие:
$$\int_{t \to \mathop t\nolimits_{раз} }^{}\mathop W\nolimits_i (t,\mathop k\nolimits_{dem} ,\mathop a\nolimits^{pr} ,m) \cdot N(t)dt \le \varepsilon > 0\;\;\;\;(18)$$где $\varepsilon$ - бесконечно малая положительная величина; $t_{раз}$ - момент разрушения (точка разрушения). При выполнении условия (18) происходит абсолютное (тотальное) разрушение. Можно ввести понятие относительного разрушения на определенном уровне:
$$\int_{t \to \mathop t\nolimits_{раз} }^{}\mathop W\nolimits_i (t,\mathop k\nolimits_{dem} ,\mathop a\nolimits^{pr} ,m) \cdot N(t)dt=\varepsilon^0 > 0\;\;\;\;(19)$$где $\varepsilon^0$ - уровень разрушения, задаваемый, например, как определенный процент $N(0)$; $t_{раз}^0$ - момент времени, в который выполняется условие (19) . Теперь перейдем к модельным экспериментам для этого решаем дифференциальное уравнение (17) при разных комбинациях параметров для всех трех сценариев, добиваясь сформулированных выше условий.
Используя «вариант а» включения $k_{dem}$ в модель, мы получили решение дифференциального уравнения (17) на диапазоне $t∈[0,43]$ при начальных условия $N(0)=221$ млн.га. Это площадь сельхозугодий России в 2009 году (Росстат, статистический сборник «Россия и страны мира 2010»). Решения в основном соотвествуют сформулированным выше условиям и получены при следующих значениях параметров $k_dem=0.9$, $a^{pr}=0.5$,$=1.45$
#=====Модель трансформации природного капитала АПК===============
#1. Позитивное воздейcтвие человека отсутствует
#2. Рассматривается первая "стратегия природы"- случай a)
#3. Рассматриваются все три функции относительной экологической нагрузки Dотн(t)
options(warn=-1)
library(deSolve) # решение дифф. уравнений с начальными условиями
library(openxlsx)
library(writexl)
library(ggthemes)
library(ggplot2)
#ВХОДНЫЕ ПАРАМЕТРЫ =============================================
at1<<-0
at2<<-0
at3<<-0
tn <<-43 #Конец сценария
a_pr <<-0.5
m <<-1.45
kdem <<-0.9
N0 <<-221
t0=seq (0,tn,1.0)
y0 <-N0
f_N2<- function (t,y, parms ){
D_otn12 <-0.2289*exp(0.0344*t)
at2=-(a_pr*D_otn12**m)+1
w <-(1 - (kdem/(at2)))
dy.dt <-y [1]*w
return (list (dy.dt ))}
out <- ode (y=y0,t=t0,func =f_N2,parms = NULL )
outN2 <- data.frame(out)
#print(outN2)
f_N3<- function (t,y, parms ){
D_otn13 <-(0.6329*exp(0.075*t))/(2.8+1.165*(-1+exp(0.075*t)))
at3=-(a_pr*D_otn13**m)+1
w <-(1 - (kdem/(at3)))
dy.dt <-y [1]*w
return (list (dy.dt ))}
out <- ode (y=y0,t=t0,func =f_N3,parms = NULL )
outN3 <- data.frame(out)
#print(outN3)
f_N1<- function (t,y, parms ){
D_otn11 <-0.229*exp(0.06*t)
at1=-(a_pr*D_otn11**m)+1
w <-(1 - (kdem/(at1)))
dy.dt <-y [1]*w
return (list (dy.dt ))}
t1=seq (0,25,1.0)
out <- ode (y=y0,t=t1,func =f_N1,parms = NULL )
outN1 <- data.frame(out)
#print(outN1)
options(repr.plot.width =14, repr.plot.height =10)
g102 <-ggplot(data=outN3, aes(x=time,y=X1))
g102 <-g102 + geom_line(data=outN3,aes(x=time,y=X1),
colour="green",linetype=1,size=1.2)
g102 <-g102 + geom_line(data=outN2,aes(x=time,y=X1),
colour="blue",linetype=1,size=1.2)
g102 <-g102 + geom_line(data=outN1,aes(x=time,y=X1),
colour="red",linetype=1,size=1.2)
g102 <-g102 +geom_text(aes(label="СЦЕНАРИЙ 1", 15,N0*0.5),
size=6,colour="red")
g102 <-g102 +geom_text(aes(label="СЦЕНАРИЙ 2", 23,N0*0.5),
size=6,colour="blue")
g102 <-g102 +geom_text(aes(label="СЦЕНАРИЙ 3", 34,N0*0.5),
size=6,colour="green")
g102 <-g102 +geom_text(aes(label="Nо=221 млн.га", 6,220),
size=6,colour="black")
g102 <-g102 +geom_text(aes(label="ЗОНА РАЗРУШЕНИЯ 10% от N0", 10,12.0),
size=4,colour="red")
g102 <-g102 + geom_hline(yintercept =N0*0.1,color="red",linetype=5, size=1.5)
g102 <-g102 +labs(x = "t",y="Природный потенциал млн.га",
title = "Сценарии развития N в отсутствии инвестиций\n (первая стратегия природы-случай а)")
#g102 <-g102 + theme_stata() + scale_colour_stata()
#g102 <-g102 +theme_economist() + scale_colour_economist()
g102 <-g102+theme_gdocs()
#x11()
print(g102)
#запись данных в excel=======================================================
r1=tn+1
r2=nrow(outN1)
dr=r1-r2
for(i in 1:dr){
outN1<- rbind(outN1,c(0,0))
}
NN.out =cbind(outN3,outN2$X1,outN1$X1)
NN.out<-data.frame(t=NN.out$time,N3=NN.out$X1,N2=outN2$X1,N1=outN1$X1)
#NN.out
Nbx <-read.xlsx("D:/#Var_prog/Var_programs/cap_N.xlsx", sheet = 1)
Nb.out=data.frame(t=Nbx$t,N1=Nbx$N1,N2=Nbx$N2,N3=Nbx$N3)
#dataset_names <- list('capN1' = Nb.out, 'capN2' = NN.out)
#openxlsx::write.xlsx(dataset_names, file = 'D:/#Var_prog/Var_programs/cap_N.xlsx')
График показывает следующие результаты. При первой «стратегии» природы (вариант а) можно ожидать, что к t=22 по сценарию 1 потери природного капитала АСЭЭС России достигнут болеее 90 процентов и система войдет в зону 10 процентного разрушения.По второму сценарию система войдет в зону разрушения природного капитала к t=34, а по третьему на конец сценареого периода t=43 уровень природного капитала составит 32.15 млн. га.При этом даже по первому сценарию первые 7 лет наблюдается рост до 250.8 млн.га, а затем начинается ускоряющися спад. По второму и третьему сценариям рост идет до 285.0 млн.га при t=11 и 280.44 млн.га при t=10 соотвественно.То есть до t=17 второй сценарий показывает себя немного лучше, чем третий.Дальше второй сценарий показывает гораздо более быстрое снижение природного капитала системы по сравнению с третьем сценарием.Полное разрушение природного капитала по сценарию 1 наступает примерно в точке t=26 остаточный уровень меньше 0.01 процента от N(0), если анализировать ситуацию дальше, то в точке приблизительно t=29 система попадает в отрицательную зону и переходит в затухающий колебательный режим.
Важно оценить, как изменяется решение от значения k_dem . На рисунке ниже показаны решения дифференциального уравнения (17)при значениях $k_{dem}$, когда позитивная роль природы отсутствует $k_{dem}=1.0$ до почти предельного позитивного значения $k_{dem}=0.1$ $(a^{pr}=0.5,m=1.45)$. Из рисунка видно, что дифференциация решений в зоне $0.5<k_{dem}≤1.0$ относительно небольшая. Затем при $k_{dem}<0.5$ разброс решений увеличивается и происходит изменения вида функции $N(t)$ от вогнутой вниз к вогнутой вверх. Это означает, что темпы (скорость) падения природного капитала меняет характер от снижающихся во времени до возрастающих.

Рассмотрим поведение $N(t)$ при третьем сценарии. Здесь картина очень интересная.При значении $k_{dem}∈[0.80,0.81]$ имеется точка бифуркации. В этой точке происходит смена тенденций в изменении N(t) на интервале $t∈[0,43]$ от падения к стабилизации, а затем к росту. При $k_{dem}=0.81$ $\;\;N(t)\;\;$стабилизируется т.е. природа полностью компенсирует отрицательное антропогенное воздействие, оказываемое по стратегии 3. При $k_{dem}=0.78$ наступает режим постоянного почти линейного роста, а дальше при падении $k_{dem}$ система срывается в положительный коллапс. При $k_{dem}= 0.915$ система переходит в режим первоначального роста, затем падения природного капитала до уровня разрушения.При $k_{dem} = 0.94$ система входит по сценарию 3 в режим постоянного падения и смещения точки входа в зону разрушения влево. При $k_{dem}$ = 1.0 (восстанавливающая сила природы отсутствует) система будет разрушена по всем трем сценариям при $t<=21$. Таким образом, данная «стратегия» природы характеризуется высоким уровнем чувствительности к объему экологической нагрузки и величине и направлению скорости ее изменения Априорные соображения, положенные в основу сценариев 2 и 3, дают возможность предполагать, что наиболее вероятной областью значений коэффициента демпфирования, при первой стратегии природы, является интервал (0.9,0.95].
Перейдем к решению дифференциального уравнения (17) при второй стратегии природы (случай б включения $k_{dem}$ в модель). Решения для следующего набора параметров $k_{dem}=0.3$ ,$a^{pr}=0.5,m=1.45$ представлены графически (см.ниже).
#=====Модель трансформации природного капитала АПК===============
#1. Позитивное воздейcтвие человека отсутствует
#2. Рассматривается вторая "стратегия природы"- случай б)
#3. Рассматриваются все три функции относительной экологической нагрузки Dотн(t)
options(warn=-1)
library(deSolve) # решение дифф. уравнений с начальными условиями
library(openxlsx)
library(writexl)
library(ggthemes)
library(ggplot2)
#ВХОДНЫЕ ПАРАМЕТРЫ =============================================
tn <<-43 #Конец сценария
a_pr <<-0.5
m <<-1.45
kdem <<-0.3
N0 <<-221
t0=seq (0,tn,1.0)
y0 <-N0
#--------решение дифференциальных уравнений 17----------------------------
#---АПОКАЛИПСИС i=1 ------------------------------------------------------
f_N1<- function (t,y, parms ){
D_otn11 <-(0.229*exp(0.06*t))
w <-(1 - (1/(1-a_pr*kdem*D_otn11**m)))
dy.dt <-y [1]*w
return (list (dy.dt ))}
out1 <- ode (y=y0,t=t0,func =f_N1,parms = NULL )
outN1 <- data.frame(out1)
#print(outN1)
#--------------ПЛОХОЙ СЦЕНАРИЙ i=2--------------------------------------------
f_N2<- function (t,y, parms ){
D_otn12 <-(0.2289*exp(0.0344*t))
w <-(1 - (1/(1-a_pr*kdem*D_otn12**m)))
dy.dt <-y [1]*w
return (list (dy.dt ))}
out2 <- ode (y=y0,t=t0,func =f_N2,parms = NULL )
outN2 <- data.frame(out2)
#print(outN2)
#--------------ХОРОШИЙ СЦЕНАРИЙ i=3--------------------------------------------
f_N3<- function (t,y, parms ){
D_otn13 <-((0.6329*exp(0.075*t))/(2.8+1.165*(-1+exp(0.075*t))))
w <-(1 - (1/(1-a_pr*kdem*D_otn13**m)))
dy.dt <-y [1]*w
return (list (dy.dt ))}
out3 <- ode (y=y0,t=t0,func =f_N3,parms = NULL )
outN3 <- data.frame(out3)
#print(outN3)
N.out <-data.frame(t=outN1[,1],N1=outN1[,2],N2=outN2[,2],N3=outN3[,2])
#print(N.out)
Nb=N.out
#===================ПОСТРОЕНИЕ ГРАФИКОВ===========================================
otrez1 <-data.frame(X= c(5,5,5),XEND =c(10,10,10),
Y=c(50,60,70),YEND=c(50,60,70))
options(repr.plot.width =14, repr.plot.height =10)
g102 <-ggplot(data=N.out, aes(x=t,y=N1))
g102 <-g102 + geom_line(data=N.out,aes(x=t,y=N1),
colour="red",linetype=1,size=1.2)
g102 <-g102 + geom_line(data=N.out,aes(x=t,y=N2),
colour="blue",linetype=1,size=1.2)
g102 <-g102 + geom_line(data=N.out,aes(x=t,y=N3),
colour="green",linetype=1,size=1.2)
g102 <-g102 +geom_text(aes(label="СЦЕНАРИЙ 1", 17,50),
size=4,colour="black")
g102 <-g102 +geom_text(aes(label="СЦЕНАРИЙ 2", 17,60),
size=4,colour="black")
g102 <-g102 +geom_text(aes(label="СЦЕНАРИЙ 3", 17,70),
size=4,colour="black")
g102 <-g102 +geom_text(aes(label="Nо=221 млн.га", 8,220),
size=6,colour="black")
g102 <-g102 +geom_text(aes(label="ЗОНА РАЗРУШЕНИЯ 10% от N0", 10,12.0),
size=4,colour="red")
g102 <-g102 + geom_hline(yintercept =N0*0.1,color="red",linetype=5, size=1.5)
g102 <-g102 +labs(x = "t",y="Природный потенциал млн.га",
title = "Сценарии развития N в отсутствии инвестиций\n (вторая стратегия природы-случай б)")
g102 <-g102+geom_segment(data=otrez1,aes(x=X,xend=XEND,y=Y,yend=YEND),
colour=c("red","blue","green"),linetype =c(1,1,1),size=c(1.5,1.5,1.5))
#g102 <-g102 + theme_stata() + scale_colour_stata()
#g102 <-g102 +theme_economist() + scale_colour_economist()
#g102 <-g102 +theme_wsj()
#g102 <-g102 +theme_tufte()
#g102 <-g102+theme_solarized()
g102 <-g102+theme_gdocs()
#g102 <-g102+ theme_void()
#g102 <-g102+ theme_classic()
#g102 <-g102+ t theme_gray()
#x11()
print(g102)
#=============запись полученных данных в EXCEL=============================================
#write.xlsx(N.out,'D:/#Var_prog/Var_programs/cap_N.xlsx',colNames = TRUE)
#write.xlsx(N.out,
#file='D:/#Var_prog/Var_programs/cap_N.xlsx',
#sheetName='capN1', append=TRUE)
#dataset_names <- list('capN1' = N.out, 'capN2' = N.out)
#openxlsx::write.xlsx(dataset_names, file = 'D:/#Var_prog/Var_programs/cap_N.xlsx')
Решения выглядят несколько иначе, чем в первом случае. Природа более слабо реагирует на незначительную негативную экологическую нагрузку (интервал [0,20]), причем тем слабее, чем ниже скорость ее возрастания. В силу этого сценарий 3 демонстрирует совершенно другой характер $N(t)$, чем в первом случае. Он значительно приближается к сценариям 1 и 2. При увеличении экологической нагрузки позитивная роль природы возрастает в результате точки разрушения в сценариях 1 и 2 немного смещаются к концу интервала [0,43], но остаются в зоне априорных соображений.
Посмотрим, как влияет значение значения $k_{dem}$ на решение при этой стратегии природы.Результаты решения дифференциального уравнения (17) при разных значениях $k_{dem}$ для сценариев 1, 2 и 3 позволяют сделать следующие выводы. При небольших значениях $k_{dem}$ 0.05, 0.1,0,2 сценарии 2,3 демонстрируют высокую стабильность (спад N(t) порядка 20 процентов от начального уровня) и незначительные различия для t от 0 до 25-30.Сценарий 3 начинает выигравать у сценария 2 при t>25-30 и тем больше , чем больше $k_{dem}$.При $k_{dem}=0.25$ только сценарий 3 не достигает при $t=43$ зоны разрушения.При $k_{dem}=0.38$ востановительные силы природы не в состоянии компенсировать негативное влияние человека и в $t=43$ сценарий 3 демонстрирует переход системы в зону разрушения. При значениях $k_{dem}>0.5$ начинается полное разрушение системы, сначало для сценария 1(здесь система уже полностью разрушалась $N(t)=0$ приблизительно при $t=30-33$) затем для сценария 2, затем почти полностью разрушается система и при сценарии 3.
Если сравнить решения для одноименных сценариев при разных стратегиях природы при разных $k_{dem}$ , то можно констатировать:
2.при первой "стратегии" поведение $N(t)$ сценарии 2,3 при разных $k_{dem}$ отличаются друг от друга кардинально
Природа при первой «стратегии» реагирует на негативноеное поведение человека особенно при малых экологических нагрузках на порядок сильнее в положительном направлении, чем при второй "стратегии".
Таким образом, вторая «стратегия природы» отличается значительной робастностью к изменению сценария поведения человека. Включение в модель той или иной «стратегии природы» определяется дополнительной априорной информацией о эколого-экономических закономерностях взаимодействия человека и природы. Кроме того, видимо, целесообразно учитывать в модели позитивное влияние природы не как константу, а как некоторую функцию $k_{dem (t)}$.
До сих пор мы рассматривали негативное влияние деятельности человека на природный капитал АСЭЭС и то, как природа может противодействовать этому. Перейдем к вопросам описаний позитивной природо-восстанавливающей деятельности человека в аграрной сфере
Рассмотрим следующее дифференциальное уравнение:
$$\frac{{dN}}{{dt}} = \Lambda (\mathop \gamma \nolimits^{Nknow} (t),\mathop I\nolimits^N (t))\;\;\;\;(20)$$То есть, прирост природного капитала, есть от растущая функция от вложенных в его развитие инвестиций и их инновационности. При этом считаем, что прирост определяется ростом продуктивности 1 га земли, тогда аналогично тому как мы описывали падение природного капитала при росте продуктивности, имеем:
$$\Delta \mathop N\nolimits_t^{S + } = \mathop N\nolimits_t \cdot (1 - \frac{1}{{\mathop b\nolimits_t }}) \ge 0\;\;\;\;(21)$$Величина $b_t$, определяющая продуктивность земли в любой момент времени, зависит от величины инвестиций на 1 га земли и уровня их инновационности. Коэффициент роста продуктивности земли $b_t$ из естественных соображений представляет собой монотонно возрастающую функцию произведения $\frac{{\mathop I\nolimits_t^N }}{{\mathop N\nolimits_0 }} \cdot \mathop \gamma \nolimits_t^{Nknow} (\mathop {inv}\nolimits_t^N = \frac{{\mathop I\nolimits_t^N }}{{\mathop N\nolimits_0 }})\;\;\;$ область значений которой интервал $[1.0,b^{br}]$.При этом в отсутствии инвестиций он должен принимать значение равное 1 , а значение $b^{pr}$, определяющее максимальный темп прироста природного капитала АСЭЭС, при предельно возможном объеме инвестиций и уровне их инновационности. Величина $\gamma_t^{Nknow}$ увеличивает, или уменьшают объем инвестиций, что либо увеличивает темп роста продуктивности земли либо снижает. Можно предложить следующие модели роста продуктивности земли от уровня инвестиций:
$$\mathop b\nolimits_t = \frac{{\mathop b\nolimits^{br} \cdot \mathop e\nolimits^{0.1 \cdot inv \cdot \mathop \gamma \nolimits_t^{Nknow} } }}{{\mathop b\nolimits^{br} + 1.0 \cdot (\mathop e\nolimits^{0.1 \cdot inv \cdot \mathop \gamma \nolimits_t^{Nknow} } - 1)}}\;\;\;(22)$$import warnings
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import rand
warnings.simplefilter('ignore')
b_br=[1.05,1.10,1.15,1.20]
inv_gam = np.arange(1, 10, 1)
Colbbr=['g','c','y','coral']
L=[r'$b^{br}=1.05$',r'$b^{br}=1.10$',r'$b^{br}=1.15$',r'$b^{br}=1.20$']
plt.figure(figsize=(12, 8))
for i in range(4):
bt=(b_br[i]*np.exp(0.1*inv_gam))/(b_br[i]+1.0*(np.exp(0.1*inv_gam))-1)
plt.plot(inv_gam ,bt,color=Colbbr[i], linewidth=2.0)
plt.xlabel(r'$inv^N*\gamma_t^{Nknow}$тыс. р',fontsize=18)
plt.ylabel(r'$b_t$',fontsize=18)
plt.title(r'$b_t(inv_t^N*\gamma^{Nknow})$', fontsize=28)
plt.legend(L, fontsize=18)
plt.grid(True)
plt.show()
Зависимость $b_t$ от $\frac{{\mathop I\nolimits_t^N }}{{\mathop N\nolimits_0 }}$ и \gamma_t^{Nknow} графически может быть представлена следующим образом:
#ГРАФИК КОЭФФИЦИЕНТА ТЕМПОВ РОСТА ПРОДУКТИВНОСТИ ЗЕМЛИ
#В ПРОСТРАНСТВЕ ИНВЕСТИЦИЙ И ИННОВАЦИЙ
import warnings
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
import ipywidgets as widgets
warnings.simplefilter('ignore')
#def theta (t):
fig = plt.figure(figsize=(12, 10))
ax = fig.gca(projection='3d')
# Make data.
b_br=1.15
inv=np.arange(1, 10, 0.5)
gamma=np.arange(0.1, 2.0, 0.1)
inv,gamma = np.meshgrid(inv, gamma)
bt=(b_br*np.exp(0.1*inv*gamma))/(b_br+1.0*(np.exp(0.1*inv*gamma))-1)
# Plot the surface.
surf = ax.plot_surface(inv, gamma, bt, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
# Customize the z axis.
ax.set_zlim(1.0, 1.12)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
ax.set_title(r'Коэффицент темпов роста продуктивности земли-$b_t$',fontsize=20)
ax.set_xlabel(r'$inv_t^N$ тыс.руб./га',fontsize=20)
ax.set_ylabel(r'$\gamma_t^{Nknow}$',fontsize=22)
ax.set_zlabel('o.e.',fontsize=22)
# Add a color bar which maps values to colors.
fig.colorbar(surf, shrink=0.5, aspect=6)
ax.dist = 6
#ax.view_init(elev=30,azim=t)
ax.view_init(elev=30,azim=195)
#plt.draw()
#plt.pause(.0001)
plt.show()
#widgets.interact(theta , t= widgets.Play(min=0,max =360));
Уровень инновационности инвестиций, можно представить как монотонно возрастающую функцию времени, область значений которой представляет собой интервал $[\gamma_{низ}^{Nknow},\gamma_{верх}^{Nknow}]$, где $\gamma_{низ}^{Nknow}<1$ предельно допустимый низший уровень инновационности инвестиций, а $\gamma_{верх}^{Nknow}≥1$ предельно возможный на интервале исследования высший уровень инновационности. Как и для других факторов развития при $\gamma_t^{Nknow=1}$ имеет место общественно необходимый уровень (средний уровень развитых экономик). Для демонстрации возможностей предлагаемого подхода возьмем следующую функцию:
$$\mathop \gamma \nolimits_t^{Nknow} = \frac{{0.4 \cdot \mathop e\nolimits^{0.246 \cdot t} }}{{1.0 + 0.4 \cdot (\mathop e\nolimits^{0.246 \cdot t} - 1)}}\;\;\;\;(23)$$import warnings
import matplotlib.pyplot as plt
import numpy as np
warnings.simplefilter('ignore')
t = np.arange(0, 43, 1)
plt.figure(figsize=(12, 8))
gamt=0.4*np.exp(0.246*t)/(1.0+0.4*(np.exp(0.246*t)-1))
plt.plot(t ,gamt,color='black', linewidth=2.0)
plt.xlabel(r'$t$ годы',fontsize=18)
plt.ylabel(r'$\gamma_t^{Nknow}$',fontsize=18)
plt.title(r'$\gamma_t^{Nknow})$ инновационность инвестиций', fontsize=28)
plt.grid(True)
plt.show()
Экономический смысл графика следующий. На начало сценарного времени инновационность инвестиций в природный капитал АСЭЭС России составляет порядка 40 процентов от общественно необходимого уровня (уровень развитых стран) предполагается, что приблизительно через 20 лет Россия достигнет 100 процентного уровня и будет поддерживать его до конца сценарного срока. Важнейшей составляющей является модель инвестиций. Предполагается, что инвестиции на 1 га земли в повышение ее продуктивности есть монотонно возрастающая функция времени с областью значений $[inv_{нач}^N,inv_{кон}^N]$. Для демонстрации возможностей подхода принимаем $inv_{нач}^N=1$ т.р/га , а $inv_{кон}^N$=10 т.р/га, стратегию роста инвестиций в сценарный период представим в виде S-образной функции:
$$\mathop {inv}\nolimits_t^N = \frac{{10 \cdot \mathop e\nolimits^{0.12 \cdot t} }}{{10 + 1.0 \cdot (\mathop e\nolimits^{0.12 \cdot t} - 1)}}\;\;\;(24)$$import warnings
import matplotlib.pyplot as plt
import numpy as np
warnings.simplefilter('ignore')
t = np.arange(0, 43, 1)
plt.figure(figsize=(12, 8))
invt=10*np.exp(0.12*t)/(10+1.0*(np.exp(0.12*t)-1))
plt.plot(t ,invt,color='black', linewidth=2.0)
plt.xlabel(r'$t$ годы',fontsize=18)
plt.ylabel(r'$inv_t^{N}$тыс.руб./га',fontsize=18)
plt.title(r'$inv_t^{N})$ инвестиции в га. земли', fontsize=28)
plt.grid(True)
plt.show()
Таким образом, располагая построенными выше соотношениями, после несложных преобразований мы приходим к следующему дифференциально-му уравнению, описывающему динамику роста природного капитала АСЭЭС:
$$\frac{{dN}}{{dt}} = V(t)N(t)\;\;\;(25)$$Где $V(t)$ -функция, в которой интегрированы введенные выше соотношения. Решение этого уравнения графически представлено ниже. Как видно из рисунка, в отсутствие негативного воздействия на природный капитал, он при заданной стратегии роста инвестиций и их инновационности может вырасти за 40 лет почти в 20 раз. Естественно в условиях значительной антропогенной нагрузки природ-ный капитал АСЭЭС так изменяться во времени не будет. Общая динамика изменения природного капитала для разных сценариев и «стратегий» при-роды может быть получена как результат решения следующих дифференци-альных уравнений:
$$\frac{{dN}}{{dt}} = V(t) \cdot N(t) + \mathop W\nolimits_i^a (t) \cdot N(t),{\rm }\frac{{dN}}{{dt}} = V(t) \cdot N(t) + \mathop W\nolimits_i^b (t) \cdot N(t),{\rm }i = 1,3{\rm }\;\;\;\;(26)$$Графически решения этих уравнений представлены ниже.
#=====Модель трансформации природного капитала АПК===============
#1. Позитивное воздейcтвие человека учитывается в модели
#2. Рассматривается первая "стратегия природы"- случай a)
#3. Рассматриваются все три функции относительной экологической нагрузки Dотн(t)
options(warn=-1)
library(deSolve) # решение дифф. уравнений с начальными условиями
library(openxlsx)
library(writexl)
library(ggthemes)
library(ggplot2)
#ВХОДНЫЕ ПАРАМЕТРЫ =============================================
at1<<-0
at2<<-0
at3<<-0
#==============================
tn <<-43 #Конец сценария
a_pr <<-0.5
m <<-1.45
kdem <<-0.90
N0 <<-221
t0=seq (0,tn,1.0)
y0 <-N0
#==============================
b.0 <<-1.0
b.n <<-1.14
gamma.0 <<-0.5
gamma.n <<-1.4
inv.0 <<- 5.0
inv.n <<- 150.0
#================================
f_N2<- function (t,y, parms ){
D_otn12 <-0.2289*exp(0.0344*t)
at2=-(a_pr*D_otn12**m)+1
w <-(1 - (kdem/(at2)))
dy.dt <-y [1]*w
return (list (dy.dt ))}
out <- ode (y=y0,t=t0,func =f_N2,parms = NULL )
outN2 <- data.frame(out)
#print(outN2)
f_N3<- function (t,y, parms ){
D_otn13 <-(0.6329*exp(0.075*t))/(2.8+1.165*(-1+exp(0.075*t)))
at3=-(a_pr*D_otn13**m)+1
w <-(1 - (kdem/(at3)))
dy.dt <-y [1]*w
return (list (dy.dt ))}
out <- ode (y=y0,t=t0,func =f_N3,parms = NULL )
outN3 <- data.frame(out)
#print(outN3)
f_N1<- function (t,y, parms ){
D_otn11 <-0.229*exp(0.06*t)
at1=-(a_pr*D_otn11**m)+1
w <-(1 - (kdem/(at1)))
dy.dt <-y [1]*w
return (list (dy.dt ))}
t1=seq (0,25,1.0)
out <- ode (y=y0,t=t1,func =f_N1,parms = NULL )
outN1 <- data.frame(out)
#print(outN1)
options(repr.plot.width =14, repr.plot.height =10)
g102 <-ggplot(data=outN3, aes(x=time,y=X1))
g102 <-g102 + geom_line(data=outN3,aes(x=time,y=X1),
colour="green",linetype=1,size=1.2)
g102 <-g102 + geom_line(data=outN2,aes(x=time,y=X1),
colour="blue",linetype=1,size=1.2)
g102 <-g102 + geom_line(data=outN1,aes(x=time,y=X1),
colour="red",linetype=1,size=1.2)
g102 <-g102 +geom_text(aes(label="СЦЕНАРИЙ 1", 16,50),
size=4,colour="black")
g102 <-g102 +geom_text(aes(label="СЦЕНАРИЙ 2", 33,32),
size=4,colour="black")
g102 <-g102 +geom_text(aes(label="СЦЕНАРИЙ 3", 40,60),
size=4,colour="black")
g102 <-g102 +geom_text(aes(label="Nо=221 млн.га", 8,220),
size=6,colour="black")
g102 <-g102 +geom_text(aes(label="ЗОНА РАЗРУШЕНИЯ 10% от N0", 10,12.0),
size=4,colour="red")
g102 <-g102 + geom_hline(yintercept =N0*0.1,color="red",linetype=5, size=1.5)
g102 <-g102 +labs(x = "t",y="Природный потенциал млн.га",title = "Сценарии развития N в отсутствии инвестиций")
g102 <-g102+theme_gdocs()
#g102 <-g102 + theme_stata() + scale_colour_stata()
#g102 <-g102 +theme_economist() + scale_colour_economist()
#x11()
#print(g102)
#----решение дифференциального уравнения 2 - позитивная деятельность человека--
f_V<- function (t,y, parms ){
GamV <-(gamma.0*exp(0.246*t))/(gamma.n+gamma.0 *(exp(0.246*t)-1))
InvV <-(inv.n*exp(0.12*t))/(inv.n + inv.0*(exp(0.12*t)-1))
giV <-InvV*GamV
bV <-(b.n*exp(0.1*giV))/(b.n+b.0*(exp(0.1*giV)-1))
V <-(1 - (1/bV))
dy.dt <-y [1]*V
return (list (dy.dt ))}
tn <<-43 #Конец сценария
t0=seq (0,tn,1.0)
y0 <-N0
out <- ode (y=y0,t=t0,func =f_V,parms = NULL )
outV <- data.frame(out)
#print(outV)
#x11()
#plot (outV[,1] ,outV[,2], type ="l", xlab ="t", ylab ="Inv",
#col ="blue", lwd =4)
g103 <-ggplot(data=outV, aes(x=outV[,1],y=outV[,2]))
g103 <-g103 + geom_line(data=outV,aes(x=outV[,1],y=outV[,2]),
colour="green",linetype=1,size=1.2)
g103 <-g103 +labs(x = "t годы",y="Природный потенциал млн.га",
title = "Сценарии N при инветициях и отсутвии экологической нагрузки")
g103 <-g103 +geom_text(aes(label="Nо=221 млн.га", 2,N0+100),
size=6,colour="black")
#g103 <-g103 + theme_stata() + scale_colour_stata()
#g103 <-g103 +theme_economist() + scale_colour_economist()
g103 <-g103+theme_gdocs()
#x11()
print(g103)
#--------решение объединенного дифференциального уравнения ----
#---АПОКАЛИПСИС----------------------------------------------------------
f_NV1<- function (t,y, parms ){
D_otn11 <-0.229*exp(0.06*t)
w <-1 - (kdem/(1-a_pr*D_otn11**m))
GamV <-(gamma.0*exp(0.246*t))/(gamma.n+gamma.0 *(exp(0.246*t)-1))
InvV <-(inv.n*exp(0.12*t))/(inv.n + inv.0*(exp(0.12*t)-1))
giV <-InvV*GamV
bV <-(b.n*exp(0.1*giV))/(b.n+b.0*(exp(0.1*giV)-1))
V <-(1 - (1/bV))
dy.dt <-y [1]*w + y [1]*V
return (list (dy.dt ))}
t1=seq (0,25,1.0)
out <- ode (y=y0,t=t1,func =f_NV1,parms = NULL )
outNV1 <- data.frame(out)
#print(outNV1)
#--------------ПЛОХОЙ СЦЕНАРИЙ--------------------------------------------
f_NV2<- function (t,y, parms ){
D_otn12 <-0.2289*exp(0.0344*t)
at2<--(a_pr*D_otn12**m)+1
w <-(1 - (kdem/(at2)))
GamV <-(gamma.0*exp(0.246*t))/(gamma.n+gamma.0 *(exp(0.246*t)-1))
InvV <-(inv.n*exp(0.12*t))/(inv.n + inv.0*(exp(0.12*t)-1))
giV <-InvV*GamV
bV <-(b.n*exp(0.1*giV))/(b.n+b.0*(exp(0.1*giV)-1))
V <-(1 - (1/bV))
dy.dt <-y [1]*w + y [1]*V
return (list (dy.dt ))}
out <- ode(y=y0,t=t0,func =f_NV2,parms = NULL )
outNV2 <- data.frame(out)
#head(out)
#x11()
#plot (outNV2[,1] ,outNV2[,2], type ="l", xlab ="t", ylab ="N2",
#col ="red", lwd =4)
#--------------ХОРОШИЙ СЦЕНАРИЙ--------------------------------------------
f_NV3<- function (t,y, parms ){
D_otn13 <-(0.6329*exp(0.075*t))/(2.8+1.165*(-1+exp(0.075*t)))
at3=-(a_pr*D_otn13**m)+1
w <-(1 - (kdem/(at3)))
GamV <-(gamma.0*exp(0.246*t))/(gamma.n+gamma.0 *(exp(0.246*t)-1))
InvV <-(inv.n*exp(0.12*t))/(inv.n + inv.0*(exp(0.12*t)-1))
giV <-InvV*GamV
bV <-(b.n*exp(0.1*giV))/(b.n+b.0*(exp(0.1*giV)-1))
V <-(1 - (1/bV))
dy.dt <-y [1]*w + y [1]*V
return (list (dy.dt ))}
out <- ode (y=y0,t=t0,func =f_NV3,parms = NULL )
outNV3 <- data.frame(out)
#head(out3)
#print(outNV3)
#x11()
#plot (outN3[,1] ,outN3[,2], type ="l", xlab ="t", ylab ="N3",
#col ="red", lwd =4)
#NV.out <-data.frame(t=outNV1[,1],N1=outNV1[,2],N2=outNV2[,2],N3=outNV3[,2])
#NV.out <-data.frame(t=outNV2[,1],N2=outNV2[,2],N3=outNV3[,2])
#NV.out <-data.frame(t=outNV3[,1],N3=outNV3[,2])
#print(NV.out)
g104 <-ggplot(data=outNV3 , aes(x=time,y=X1))
g104 <-g104 + geom_line(data=outNV3,aes(x=time,y=X1),
colour="green",linetype=1,size=1.2)
g104 <-g104 + geom_line(data=outNV2,aes(x=time,y=X1),
colour="blue",linetype=1,size=1.2)
g104 <-g104 + geom_line(data=outNV1,aes(x=time,y=X1),
colour="red",linetype=1,size=1.2)
g104 <-g104 +geom_text(aes(label="СЦЕНАРИЙ 1", 25,outNV1[26,2]+60),
size=4,colour="red")
g104 <-g104 +geom_text(aes(label="СЦЕНАРИЙ 2", 41,outNV2[44,2]+60),
size=4,colour="blue")
g104 <-g104 +geom_text(aes(label="СЦЕНАРИЙ 3", 40,outNV3[44,2]-50),
size=4,colour="chartreuse4")
g104 <-g104 +geom_text(aes(label="Nо=221 млн.га", 2,N0-20),
size=6,colour="black")
g104 <-g104 +geom_text(aes(label="ЗОНА РАЗРУШЕНИЯ 10% от N0", 10,8.0),
size=4,colour="red")
g104 <-g104 + geom_hline(yintercept =N0*0.1,color="red",linetype=5, size=1.5)
g104 <-g104 +labs(x = "t",y="Природный потенциал млн.га",
title = "Сценарии развития N при инвестициях\n первая стратегия природы")
#g104 <-g104 + theme_stata() + scale_colour_stata()
g104 <-g104+theme_gdocs()
#x11()
print(g104)
#запись данных в excel=======================================================
r1=tn+1
r2=nrow(outN1)
dr=r1-r2
for(i in 1:dr){
outNV1<- rbind(outNV1,c(0,0))
}
NNV.out =cbind(outNV3,outNV2$X1,outNV1$X1)
NNV.out<-data.frame(t=NNV.out$time,N3=NNV.out$X1,N2=outNV2$X1,N1=outNV1$X1)
#NNV.out
#Nbx <-read.xlsx("D:/#Var_prog/Var_programs/cap_N.xlsx", sheet = 1)
#Nb.out=data.frame(t=Nbx$t,N1=Nbx$N1,N2=Nbx$N2,N3=Nbx$N3)
#dataset_names <- list('capN1' = Nb.out, 'capN2' = NN.out)
#openxlsx::write.xlsx(dataset_names, file = 'D:/#Var_prog/Var_programs/cap_N.xlsx')
#=====Модель трансформации природного капитала АПК===============
#1. Позитивное воздейcтвие человека учитывается в модели
#2. Рассматривается вторая "стратегия природы"- случай б)
#3. Рассматриваются все три функции относительной экологической нагрузки Dотн(t)
options(warn=-1)
library(deSolve) # решение дифф. уравнений с начальными условиями
library(openxlsx)
library(writexl)
library(ggthemes)
library(ggplot2)
#ВХОДНЫЕ ПАРАМЕТРЫ ============================================
#==============================
tn <<-43 #Конец сценария
a_pr <<-0.5
m <<-1.45
kdem <<-0.45
N0 <<-221
t0=seq (0,tn,1.0)
y0 <-N0
#==============================
b.0 <<-1.0
b.n <<-1.14
gamma.0 <<-0.5
gamma.n <<-1.4
inv.0 <<- 5.0
inv.n <<- 150.0
#================================
#--------решение объединенного дифференциального уравнения ----
#---СЦЕНАРИЙ 1 АПОКАЛИПСИС-----------------------------------------------------
f_NV1<- function (t,y, parms ){
D_otn21 <-0.229*exp(0.06*t)
w21 <-1 - 1/(1-a_pr*kdem*D_otn21**m)
GamV21 <-(gamma.0*exp(0.246*t))/(gamma.n+gamma.0 *(exp(0.246*t)-1))
InvV21 <-(inv.n*exp(0.12*t))/(inv.n + inv.0*(exp(0.12*t)-1))
giV21 <-InvV21*GamV21
bV21 <-(b.n*exp(0.1*giV21))/(b.n+b.0*(exp(0.1*giV21)-1))
V21 <-(1 - (1/bV21))
dy.dt <-y [1]*w21 + y [1]*V21
return (list (dy.dt ))}
out4 <- ode (y=y0,t=t0,func =f_NV1,parms = NULL )
outNV1 <- data.frame(out4)
#x11()
#plot (outNV1[,1] ,outNV1[,2], type ="l", xlab ="t", ylab ="NV1",
#col ="blue", lwd =4)
#--------------СЦЕНАРИЙ 2 ПЛОХОЙ------------------------------------------
f_NV2<- function (t,y, parms ){
D_otn22 <-0.2289*exp(0.0344*t)
w22 <-1 - 1/(1-a_pr*kdem*D_otn22**m)
GamV22 <-(gamma.0*exp(0.246*t))/(gamma.n+gamma.0 *(exp(0.246*t)-1))
InvV22 <-(inv.n*exp(0.12*t))/(inv.n + inv.0*(exp(0.12*t)-1))
giV22 <-InvV22*GamV22
bV22 <-(b.n*exp(0.1*giV22))/(b.n+b.0*(exp(0.1*giV22)-1))
V22 <-(1 - (1/bV22))
dy.dt <-y [1]*w22 + y [1]*V22
return (list (dy.dt ))}
out5 <- ode (y=y0,t=t0,func =f_NV2,parms = NULL )
outNV2 <- data.frame(out5)
#head(out)
#x11()
#plot (outNV2[,1] ,outNV2[,2], type ="l", xlab ="t", ylab ="N2",
#col ="red", lwd =4)
#--------------СЦЕНАРИЙ 3--ХОРОШИЙ СЦЕНАРИЙ-----------------------------------
f_NV3<- function (t,y, parms ){
D_otn13 <-(0.6329*exp(0.075*t))/(2.8+1.165*(-1+exp(0.075*t)))
w <-1 - 1/(1-a_pr*kdem*D_otn13**m)
GamV <-(gamma.0*exp(0.246*t))/(gamma.n+gamma.0 *(exp(0.246*t)-1))
InvV <-(inv.n*exp(0.12*t))/(inv.n + inv.0*(exp(0.12*t)-1))
giV <-InvV*GamV
bV <-(b.n*exp(0.1*giV))/(b.n+b.0*(exp(0.1*giV)-1))
V <-(1 - (1/bV))
dy.dt <-y [1]*w + y [1]*V
return (list (dy.dt ))}
tn <<-43 #Конец сценария
t0=seq (0,tn,1.0)
y0 <-N0
out <- ode (y=y0,t=t0,func =f_NV3,parms = NULL )
outNV3 <- data.frame(out)
#head(out)
#x11()
#plot (outN3[,1] ,outN3[,2], type ="l", xlab ="t", ylab ="N3",
#col ="red", lwd =4)
NV.out <-data.frame(t=outNV1[,1],N1=outNV1[,2],N2=outNV2[,2],N3=outNV3[,2])
#print(NV.out)
options(repr.plot.width =14, repr.plot.height =10)
g104 <-ggplot(data=NV.out, aes(x=t,y=N1))
g104 <-g104 + geom_line(data=NV.out,aes(x=t,y=N1),
colour="red",linetype=1,size=1.2)
g104 <-g104 + geom_line(data=NV.out,aes(x=t,y=N2),
colour="blue",linetype=1,size=1.2)
g104 <-g104 + geom_line(data=NV.out,aes(x=t,y=N3),
colour="green",linetype=1,size=1.2)
g104 <-g104 +geom_text(aes(label="СЦЕНАРИЙ 1", 17,outNV1[29,2]+60),
size=4,colour="red")
g104 <-g104 +geom_text(aes(label="СЦЕНАРИЙ 2", 34,outNV2[44,2]+40),
size=4,colour="blue")
g104 <-g104 +geom_text(aes(label="СЦЕНАРИЙ 3", 38,outNV3[44,2]-50),
size=4,colour="chartreuse4")
g104 <-g104 +geom_text(aes(label="Nо=221 млн.га", 10,N0+10),
size=5,colour="black")
g104 <-g104 +geom_text(aes(label="ЗОНА РАЗРУШЕНИЯ 10% от N0", 10,12.0),
size=4,colour="red")
g104 <-g104 + geom_hline(yintercept =N0*0.1,color="red",linetype=5, size=1.5)
g104 <-g104 +labs(x = "t",y="Природный потенциал млн.га",
title = "Сценарии развития N при инвестициях\n вторая стратегия природы")
#g104 <-g104 + theme_stata() + scale_colour_stata()
#g104 <-g104 +theme_economist() + scale_colour_economist()
g104 <-g104+theme_gdocs()
print(g104)
Графики наглядно демонстрируют, что негативные сценарии (1,2) по экологической нагрузке не могут быть компенсированы только за счет инвестиций в повышение продуктивности земли при любой стратегии природы, то есть разрушение практически неизбежно.А при экологичном сценарии (3) в зависимости от стратегии природы может быть обеспечен и значительный рост за 43 года с дальнейшей стабилизацией(вариант а) или снижение до 40-50 процентов с дальнешем ростом до первоначального уровня с небольшим ростом и стабилизацией (вариант б) Полученные результаты по формализации описания динамики основного, человеческого и природного капитала показывают достаточно большие возможности предлагаемого подхода для сценарного анализа развития эко-номического потенциала АСЭЭС на больших временных горизонтах. Необходимо обобщение полученных результатов в виде достаточно строгой теории трансформации факторов развития АСЭЭС.
library(openxlsx)
library(writexl)
#dataset_names <- list('WVb_01' = dat025,'WVb_02' = dat030,'WVb_03' = dat035,'WVb_04' = dat040,'WVb_5' = dat045,
#'WVb_06' = dat050,'WVb_07' = dat055,'WVb_08' = dat060)
#openxlsx::write.xlsx(dataset_names, file = 'CapN_Data.xlsx')