主体的收益函数¶
%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
import matplotlib.pyplot as plt
from pathlib import Path
import numpy as np
import pandas as pd
社会收益函数¶
参考 [@castillarho2017a] 的研究方法，可以使用自身违反规范的次数和道格拉斯函数（Cobb-Douglas Function）作为收益函数，评估某项社会成本 $cost$ 随观测次数 $x$ 的非线性变化。
\begin{equation} s_{x} = 1 - \beta^x \end{equation}
下图展示了该函数随不同参数 $\beta$ 和变量 $x$ 取值不同而产生的成本变化。
from src.func import cobb_douglas
x = np.arange(6)
for p in np.arange(0, 1, 0.1):
y = np.vectorize(cobb_douglas)(p, x)
plt.plot(x, y, label=f"$\\beta={p:.2f}$")
# 添加注释
plt.xlabel(r"$x$")
plt.ylabel("$cost$ [0, 1]")
plt.legend()
plt.show();
在我们的模型设定中，一个具体主体主体 $v$ 的社会得分 $S_v$ 包括两部分：
- $s_m$: 违反制度而导致自身的声誉损失
- $s_n$: 由于其他主体违反制度而产生不满
两部分都可以同样使用上述柯布-道格拉斯函数计算进行计算，其中 $m$ 是主体 $v$ 的违规次数，$n$ 是主体 $v$ 对其他主体违规的不满次数。因为这两部分的社会成本是相互独立的。最后，我们将两部分的社会成本进行等权重平均：
\begin{equation} S_v = (s_{m} + s_{n}) / 2 \end{equation}
至于使用式(1)所示的柯布-道格拉斯方程计算 $s_m$ 和 $s_n$ 时的两个参数 $\beta_m$ 和 $\beta_n$，我们同样参考前人研究，基于社会价值观问卷调查的估计：
- $\beta_1$: 反映主体对于自身遵循规范的严格程度，如果过于严格，则违反制度的行为会严重降低声誉。
- $\beta_2$: 反映主体倾向于和集体决策保持一致的程度，如果过于集体主义，则会对其他人不一致的违规行为感到厌恶。
from src.func import lost_reputation
from hydra import compose, initialize
import seaborn as sns
import os
# 加载项目层面的配置
with initialize(version_base=None, config_path="../../config"):
cfg = compose(config_name="config")
os.chdir(cfg.root)
# ------ 绘制图片 --------
m = np.arange(6, dtype=int)
n = np.arange(6, dtype=int)
x, y = np.meshgrid(m, n)
beta1 = 1 - cfg.farmer.s_enforcement_cost
beta2 = cfg.farmer.s_reputation
s = lost_reputation(
cost=beta1,
reputation=beta2,
caught_times=x,
punish_times=y,
)
print(rf"Current default $\beta_1$: {beta1}, $\beta_2$: {beta2}.")
sns.heatmap(s, annot=True)
plt.xlabel("Breach caught times")
plt.ylabel("Less tolerant to non-compliance")
plt.show();
Current default $\beta_1$: 0.74, $\beta_2$: 0.55.
经济收益函数¶
主体通过灌溉生产粮食，根据粮食价格和水资源灌溉成本，取得经济收益。
\begin{equation} E_{v} = p_{c} * Y_{c, v} - V_{sw, v} * p_{sw} - V_{gw, v} * p_{gw} \end{equation}
其中 $c$ 是主体 $v$ 生产的作物类型，$p_{c}$ 为该作物的市场价格（$CNY/t$），$Y_{c, v}$ 是该主体生产此作物的总产量（单位：$t$）。$p_{sw}$ 和 $p_{gw}$ 分别为该主体所在地区的地表水、地下水价格（单位 $CNY/m^3$），$V_{sw, v}$ 和 $V_{gw, v}$ 则分别为该主体的总灌溉用水量（单位 $m^3$）。
水价计算¶
地下水价格根据则根据机井的电费和抽水量计算，不同地区根据地下水位不同差别很大，估计时一般用“以电折水”的办法进行计算。
参考相关研究文章，我们设定地下水的价格为统一的平均水平 $0.68$ 元/立方米。
综上，我们使用的各省地表/地下水价格数据如下表所示：
# 加载地表、地下水价格数据
water_prices = pd.read_csv(cfg.ds.water_prices, index_col=0)
water_prices
| surface | ground | |
|---|---|---|
| name_en | ||
| Henan | 0.4500 | 0.68 |
| Neimeng | 0.2600 | 0.68 |
| Ningxia | 0.0305 | 0.68 |
| Qinghai | 0.3700 | 0.68 |
| Shandong | 0.2400 | 0.68 |
| Shanxi | 0.2500 | 0.68 |
| Shaanxi | 0.1900 | 0.68 |
| Gansu | 0.3400 | 0.68 |
粮食生产价格数据来自于国家统计局农村社会经济调查司的中国农产品价格调查年鉴，针对三种主要作物有价格差异（单位：元/吨）。
cfg.ds.crop_prices
{'Rice': 8000, 'Wheat': 4000, 'Maize': 2330}
那么，对于一个灌溉面积为 1 的单元，其经济收益随地表/地下水变化的趋势大致为：
from functools import partial
from aquacrop_abses.payoff import economic_payoff
# ----- setup 一个例子 ----
area = 1 # 灌溉面积（公顷）
q_surface = np.arange(100) # 地表水灌溉量（毫米）
q_ground = np.arange(100) # 地下水灌溉量（毫米）
x, y = np.meshgrid(q_surface, q_ground)
province = "Henan"
prices = water_prices.loc[province]
# 计算用水成本
func = partial(economic_payoff, water_prices=prices)
costs = np.vectorize(func)(q_surface=x, q_ground=y)
# ---- 绘制用水成本 ----
_, ax = plt.subplots()
sns.heatmap(costs, ax=ax)
ax.set_xlabel("Surface water price")
ax.set_ylabel("Ground water price")
plt.show();
综合收益函数¶
最终，我们将主体 $v$ 的社会收益和经济收益分别在其所在的社会联系中进行比较，将两类收益都转化成排名，并对排名进行加权平均，得到主体 $v$ 的综合收益指标 $R_v$（0～1）：
\begin{equation} R_v = \frac{f(S_v) + f(E_v)}{2} \end{equation}
其中：
\begin{equation} f(R) = \frac{R_{\max} - R_{\min}}{R_v - R_{\min}} \end{equation}
上式中 $R_{\max} = \max{R_{v'}}$，$v' \in L_v$，$L_v$ 代表与主体 $v$ 相关联的所有主体集合。