示例: Tumurbulag Soum 模拟计算
==============================

这里以 Tumurbulag Soum 的计算为例，说明如何在 Jupyter 中进行模拟计算。

数据准备与计算
--------------

.. code:: python3

    import pysd
    import numpy as np
    import pandas as pd
    import re
    from pathlib import Path
    import matplotlib.pyplot as plt

    import warnings
    warnings.filterwarnings('ignore')

定义三种情景模拟参量。

.. code:: python3

    keys = ['Forage production','Carrying capacity','Livestock','Livestock output','Livestock output rate']

定义函数，针对多情景模拟进行计算。

.. code:: python3

    def calculate(project_path):

        historical_livestock = pd.read_excel( project_path / 'input_data.xlsx',
                                             'historical_livestock').iloc[1][2:].values.astype('float')
        for scenario in [126,245,585]:
            fn = project_path / 'models' / f'{scenario}.py'
            strings = open(fn).read()
            modified_string = re.sub(r',\n    lambda: (\d+)', r',\n    lambda: '+str(int(historical_livestock[0])), strings)
            open(fn,'w',encoding='utf-8').write(modified_string)

        x = np.log(np.arange(1,historical_livestock.shape[0]+1))
        BAUk,BAUb = np.polyfit(x,historical_livestock,1)
        x = np.arange(36)

        BAU_values = BAUk*np.log(x)+BAUb

        all_results = []
        for scenario in [126,245,585]:
            model = pysd.load( project_path / 'models' / f'{scenario}.py' )
            all_results.append(model.run()[keys])
        for key in keys:
            all_results[0][key].iloc[7] = all_results[1].iloc[7][key]
            all_results[1][key].iloc[7] = all_results[1].iloc[7][key]
            all_results[2][key].iloc[7] = all_results[1].iloc[7][key]

        all_results[0]['Livestock'].iloc[7] = historical_livestock[-1] # 126
        all_results[1]['Livestock'].iloc[7] = historical_livestock[-1] # 245
        all_results[2]['Livestock'].iloc[7] = historical_livestock[-1] # 585

        return all_results,BAU_values

2022-2050年畜牧业系统情景
-------------------------

.. code:: python3

    project_path = Path('Tumurbulag')
    all_results,  BAU_values = calculate(project_path)
    results = [i[7:] for i in all_results]

定义函数，对模拟数据进行制图。

.. code:: python3

    def show_results(results):
        colors = ['blue','red','green']
        plt.figure(figsize=(16,10))
        plt.subplot(2,2,1)
        for i in range(3):
            plt.plot(np.arange(2022,2051),results[i]['Forage production'].values,colors[i])
            plt.legend(['SSP126','SSP245','SSP585'],loc='upper right')
        plt.title('Scenarios of Forage production during 2022-2050 ')

        plt.subplot(2,2,2)
        for i in range(3):
            plt.plot(np.arange(2022,2051),results[i]['Carrying capacity'].values,colors[i])
            plt.legend(['SSP126','SSP245','SSP585'],loc='upper right')
        plt.title('Scenarios of Carrying capacity during 2022-2050 ')

        plt.subplot(2,2,3)
        for i in range(3):
            plt.plot(np.arange(2022,2051),results[i]['Livestock'].values,colors[i])
            plt.legend(['SSP126','SSP245','SSP585'],loc='upper right')
        plt.title('Scenarios of Livestok during 2022-2050 ')

        plt.subplot(2,2,4)
        for i in range(3):
            plt.plot(np.arange(2022,2051),results[i]['Livestock output'].values,colors[i])
            plt.legend(['SSP126','SSP245','SSP585'],loc='upper right')
        plt.title('Scenarios of Livestock output during 2022-2050 ')
    show_results(results)



.. image:: ./imgs_files/imgs_9_0_xac.png


定义函数计算调控参数：

.. code:: python3

    def scen(results):
        means = []
        stds = []
        MV = []
        SV = []
        for key in keys:
            mean_ = np.array([re[key].mean() for re in results])
            MV.extend(mean_)
            std_ = np.array([re[key].std() for re in results])
            SV.extend(std_)
            means.append((mean_==(mean_.max()))+(mean_==(mean_.min()))*-1)
            stds.append((std_==(std_.min()))+(std_==(std_.max()))*-1)
        MV.extend([np.nan]*3)
        SV.extend([np.nan]*3)

        index = pd.MultiIndex.from_product(
            [["Forage production/ton", "Carrying capacity/SU","Livestock/SU","Livestock output/SU","Livestock output rate/SU","Sustainability scores"],
             ["SSP126", "SSP245", "SSP585"]],
            names=["Variable", "Scenarios"]
        )

        columns = ["MV", "SV", "MV'", "SV'", "TV_j"]

        MV_ = np.concatenate([np.array(means).flatten(),np.array(means).sum(axis=0)])
        SV_ =np.concatenate([np.array(stds).flatten(),np.array(stds).sum(axis=0)])
        TV = MV_+SV_
        data = np.concatenate([[MV],[SV],[MV_],[SV_],[TV]])
        Sus_df = pd.DataFrame(data.T, index=index, columns=columns)
        return Sus_df, means, stds


进行计算，查看结果：

.. code:: python3

    Sus_df,means, stds = scen(results)
    Sus_df




定义函数，计算调控结果：

.. code:: python3

    def show_text():

        TV = np.array(means).sum(axis=0)+np.array(stds).sum(axis=0)
        if TV[1]>=2:
            SSP_index=1
        else:
            SSP_index = np.argmax(TV)
        Best_Scenario = 'SSP'+str([126,245,585][SSP_index])

        text = f"""Best Scenario is {Best_Scenario}.
        Over the next three decades, forage production, carrying capacity, livestock numbers and livestock output fluctuate and remain stable.
        Among the three scenarios (SSP126, SSP245 and SSP585),the {Best_Scenario} scenario for {project_path.name} in the lower reaches of the basin,
        was identified as the rational scenario for the highest productivity and stability.
        """
        print(text)
        return SSP_index, Best_Scenario


运行上面函数：

.. code:: python3

    SSP_index, Best_Scenario = show_text()


.. parsed-literal::

    Best Scenario is SSP126.
        Over the next three decades, forage production, carrying capacity, livestock numbers and livestock output fluctuate and remain stable.
        Among the three scenarios (SSP126, SSP245 and SSP585),the SSP126 scenario for Tumurbulag in the lower reaches of the basin,
        was identified as the rational scenario for the highest productivity and stability.



历史期和预测期之间的库存率变化
------------------------------

定义函数，计算变化 ：

.. code:: python3

    def change():
        historical_data = pd.read_excel(project_path / 'input_data.xlsx','grassland_area').iloc()[SSP_index+1][2:9]
        grass_land_area = np.nanmean(historical_data.values)

        Carrying_capacity_Scenario = Sus_df['MV']['Carrying capacity/SU'][SSP_index]
        Forage_production_Scenario = Sus_df['MV']['Forage production/ton'][SSP_index]
        Livestock_Scenario = Sus_df['MV']['Livestock/SU'][SSP_index]
        Livestock_output_Scenario = Sus_df['MV']['Livestock output/SU'][SSP_index]

        Forage_production_before = all_results[SSP_index][:8]['Forage production'].mean()
        Carrying_capacity_before = all_results[SSP_index][:8]['Carrying capacity'].mean()
        Livestock_before = all_results[SSP_index][:8]['Livestock'].mean()
        Livestock_output_before = all_results[SSP_index][:8]['Livestock output'].mean()


        data = [[project_path.name,
                 '2015-2022',
                 '%d'%(Forage_production_before/grass_land_area*1000),
                 '%d' %(Carrying_capacity_before),
                 '%.2f'%(Carrying_capacity_before/grass_land_area),
                 '%d'%(Livestock_before),'%.2f'%(Livestock_before/Carrying_capacity_before),
                 '%d'%(Livestock_output_before)],
                ['',
                 '2022-2050',
                 '%d'%(Forage_production_Scenario/grass_land_area*1000),
                 '%d' %(Carrying_capacity_Scenario),
                 '%.2f'%(Carrying_capacity_Scenario/grass_land_area),
                 '%d'%(Livestock_Scenario),
                 '%.2f'%(Livestock_Scenario/Carrying_capacity_Scenario),
                 '%d'%(Livestock_output_Scenario)]
               ]

        # Define indexes for multiple layers of columns
        columns = pd.MultiIndex.from_tuples([
            ('Soums', ''),
            ('Period', ''),
            ('Grassland', 'Forage (kg/ha)'),
            ('Grassland', 'Carrying capacity (SU)'),
            ('Grassland', 'Carrying capacity (SU/ha)'),
            ('Livestock', 'Livestock (SU)'),
            ('Livestock', 'Stocking rate'),
            ('Livestock', 'Livestock output (SU)')
        ])
        change_df = pd.DataFrame(data, columns=columns)

        num1 = (Forage_production_Scenario-Forage_production_before)*100/Forage_production_before
        num2 = Carrying_capacity_before/grass_land_area
        num3 = Carrying_capacity_Scenario/grass_land_area
        num4 = Livestock_before/Carrying_capacity_before
        num5 = Livestock_Scenario/Carrying_capacity_Scenario
        text = f'''The projected average forage yield in {project_path.name} is expected to {'increase' if num1 > 0 else'decrease'} by {'%.2f'%(abs(num1))}%,
        compared to the historical period.
        The carrying capacity is projected to {'decrease' if num2<num3 else 'increase'} from {'%.2f'%num2} SU/ha during 2015-2022 to {'%.2f'%num3} SU/ha during 2022-2050.
        The stocking rate is projected to {'decrease' if num4<num5 else 'increase'} from {'%.2f'%num4} during 2015-2022 to {'%.2f'%num5} during 2022-2050.'''
        print(text)
        return change_df


运行函数，生成结果：

.. code:: python3

    change_df = change()


.. parsed-literal::

    The projected average forage yield in Tumurbulag is expected to increase by 0.67%,
        compared to the historical period.
        The carrying capacity is projected to decrease from 1.19 SU/ha during 2015-2022 to 1.21 SU/ha during 2022-2050.
        The stocking rate is projected to decrease from 1.46 during 2015-2022 to 1.55 during 2022-2050.


.. code:: python3

    change_df




畜牧系统回归模拟
----------------

定义函数，计算根据情景模拟结果的策略建议。

.. code:: python3

    def strategy():
        BAU = BAU_values[8:].mean()
        Cp = results[SSP_index].mean()['Carrying capacity']
        stock_rate_BAU = BAU/Cp
        SSP_RCP_Scenarios = results[SSP_index].mean()['Livestock']
        Stock_rate_Scenario = SSP_RCP_Scenarios/Cp
        BAU_RCP_Reduction = BAU-SSP_RCP_Scenarios
        reduction_rate_0 = BAU_RCP_Reduction/BAU
        RCP_Cp_Reduction = SSP_RCP_Scenarios-Cp
        reduction_rate_1 = RCP_Cp_Reduction/SSP_RCP_Scenarios

        columns = "Livestock	BAU	SSP-RCP Scenarios	Carrying capacity	Reduction from BAU to rational scenarios	Reduction rate from BAU to rational scenarios	Reduction from scenarios to Carrying capacity	Reduction rate from scenarios to Carrying capacity ".split('\t')
        data = [project_path.name,BAU,SSP_RCP_Scenarios,Cp,BAU_RCP_Reduction,reduction_rate_0,RCP_Cp_Reduction,reduction_rate_1]
        strategy_df = pd.DataFrame(np.array([data]), columns=columns)
        text = f"""In the coming three decades, {project_path.name} will be heavily overloaded under the BAU scenario,
        and even under the rational SSP scenarios.
        To maintain forage-livestock equilibrium,livestock should be reduced in the two consecutive steps:
        1) First Step:reducing livestock from BAU torational future scenarios by {'%d' %BAU_RCP_Reduction} SU ({'%.2f' %(reduction_rate_0*100)}% reduction).
        2) SecondStep: reducing livestock from rational scenarios to carrying capacity {'%d' %RCP_Cp_Reduction} SU ({'%.2f' %(reduction_rate_1*100)}% reduction)."""
        print(text)

        return strategy_df


运行函数，查看结果：

.. code:: python3

    strategy_df = strategy()
    strategy_df


.. parsed-literal::

    In the coming three decades, Tumurbulag will be heavily overloaded under the BAU scenario,
        and even under the rational SSP scenarios.
        To maintain forage-livestock equilibrium,livestock should be reduced in the two consecutive steps:
        1) First Step:reducing livestock from BAU torational future scenarios by 171060 SU (26.89% reduction).
        2) SecondStep: reducing livestock from rational scenarios to carrying capacity 164164 SU (35.31% reduction).




定义函数，进行拟合，并绘图：

.. code:: python3

    def plot_fit():
        historical_livestock = pd.read_excel( project_path / 'input_data.xlsx',
                                             'historical_livestock').iloc[1][2:].values.astype('float')
        plt.plot(np.arange(2015,2023),historical_livestock,color='black')
        plt.plot(np.arange(2022,2051),results[SSP_index]['Livestock'].values,color='#0000FF')
        x,y = [np.arange(2022,2051),results[SSP_index]['Livestock'].values]
        k,b = np.polyfit(x,y,1)

        print(k,b)

        plt.plot(np.arange(2015,2051),k*np.arange(2015,2051)+b,'.',color='#4472C4',markersize=3)
        x,y = [np.arange(2022,2051),results[SSP_index]['Carrying capacity'].values]
        k,b = np.polyfit(x,y,1)
        print(k,b)

        plt.plot(x,y,color='#7F7F7F')
        plt.plot(np.arange(2015,2051),k*np.arange(2015,2051)+b,'.',color='#7F7F7F',markersize=3)
        plt.plot(np.arange(2015,2051),BAU_values,'--',c='black')
        plt.xlim(2015,2050)
        plt.legend(['Historical livestock',Best_Scenario,f'Linear({Best_Scenario})','Carrying capacity','Linear(Carrying capacity)','BAU'],
                   bbox_to_anchor=(1.05,0),
                   loc='lower left')


运行此函数：

.. code:: python3

    plot_fit()


.. parsed-literal::

    4104.998173025084 -7892787.658744299
    3737.161682109544 -7308037.384229856



.. image:: ./imgs_files/imgs_30_1_xk3.png