{ "cells": [ { "cell_type": "markdown", "id": "8a57cf99-9019-4022-9546-4f86b1d9ca88", "metadata": { "ExecuteTime": { "end_time": "2022-09-30T21:49:22.516232Z", "start_time": "2022-09-30T21:49:22.430418Z" } }, "source": [ "```{index} single: application; resource allocation\n", "```\n", "```{index} single: solver; cbc\n", "```\n", "```{index} pandas dataframe\n", "```\n", "```{index} single: AMPL; sets\n", "```\n", "```{index} stochastic optimization\n", "```\n", "```{index} two-stage problem\n", "```\n", "\n", "# Extra: The farmer's problem and its variants\n", "\n", "The [Farmer's Problem](https://www.math.uh.edu/~rohop/Spring_15/Chapter1.pdf) is a teaching example presented in the well-known textbook by John Birge and Francois Louveaux.\n", "\n", "* Birge, John R., and Francois Louveaux. Introduction to stochastic programming. Springer Science & Business Media, 2011." ] }, { "cell_type": "code", "execution_count": null, "id": "9720e9d1", "metadata": { "ExecuteTime": { "end_time": "2022-09-30T21:49:20.540447Z", "start_time": "2022-09-30T21:49:20.292334Z" } }, "outputs": [], "source": [ "# install dependencies and select solver\n", "%pip install -q amplpy numpy pandas matplotlib\n", "\n", "SOLVER = \"cbc\"\n", "\n", "from amplpy import AMPL, ampl_notebook\n", "\n", "ampl = ampl_notebook(\n", " modules=[\"cbc\"], # modules to install\n", " license_uuid=\"default\", # license to use\n", ") # instantiate AMPL object and register magics" ] }, { "cell_type": "code", "execution_count": 2, "id": "330cd1f1", "metadata": { "ExecuteTime": { "end_time": "2022-09-30T21:49:22.291926Z", "start_time": "2022-09-30T21:49:20.543823Z" } }, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "4bdf92e9-c0b8-407e-8db3-817896abeadd", "metadata": { "ExecuteTime": { "end_time": "2022-09-30T21:49:22.516232Z", "start_time": "2022-09-30T21:49:22.430418Z" } }, "source": [ "## Problem Statement\n", "\n", "In the [farmer's problem](https://www.math.uh.edu/~rohop/Spring_15/Chapter1.pdf), a European farmer must allocate 500 acres of land among three different crops (wheat, corn, and sugar beets) aiming to maximize profit. The following information is given:\n", "\n", "* Planting one acre of wheat, corn and beet costs 150, 230 and 260 euro, respectively.\n", "\n", "* The mean yields are 2.5, 3.0, and 20.0 tons per acre for wheat, corn, and sugar beets, respectively. The yields can vary up to 20% from nominal conditions depending on weather.\n", "\n", "* At least 200 tons of wheat and 240 tons of corn are needed for cattle feed. Cattle feed can be raised on the farm or purchased from a wholesaler. \n", "\n", "* The mean selling prices have been 170 and 150 euro per ton of wheat and corn, respectively, over the last decade. The purchase prices are 40% more due to wholesaler's margins and transportation costs.\n", "\n", "* Sugar beets are a profitable crop expected to sell at 36 euro per ton, but there is a quota on sugar beet production. Any amount in excess of the quota can be sold at only 10 euro per ton. The farmer's quota for next year is 6,000 tons.\n", "\n", "Create three solutions for the farmer to consider.\n", "\n", "1. How should the farmer allocate land to maximize profit for the mean weather conditions?\n", "\n", "2. The second solution should explicitly consider the impact of weather. How should the farmer allocate land to maximize expected profit if the yields could go up or down by 20% due to weather conditions? What is the profit under each scenario?\n", "\n", "3. During your interview you learned the farmer needs a minimal profit each year to stay in business. How would you allocate land use to maximize the worst case profit? \n", "\n", "4. Determine the tradeoff between risk and return by computing the mean expected profit when the minimum required profit is the worst case found in part 3, and 58,000, 56,000, 54,000, 52,000, 50,000, and 48,000 euro. Compare these solutions to part 2 by plotting the expected loss in profit. \n", "\n", "5. What would be your advice to the farmer regarding land allocation?" ] }, { "cell_type": "markdown", "id": "a5555bbd-5865-412a-8933-eceab4f29ceb", "metadata": { "ExecuteTime": { "end_time": "2022-09-30T21:49:22.516232Z", "start_time": "2022-09-30T21:49:22.430418Z" } }, "source": [ "## Data Summary\n", "\n", "| Scenario | Yield for wheat
(tons/acre)| Yield for corn
(tons/acre) | Yield for beets
(tons/acre) |\n", "| :-- | :-: | :-: | :-: |\n", "| Good weather | 3 | 3.6 | 24 |\n", "| Average weather | 2.5 | 3 | 20 |\n", "| Bad weather | 2 | 2.4 | 16 |\n", "\n", "We first consider the case in which all the prices are fixed and not weather-dependent. The following table summarizes the data.\n", "\n", "| Commodity | Sell
Price
(euro/ton) | Market
Demand
(tons) | Excess Price | Buy
Price
(euro/ton) | Cattle Feed
Required
(tons) | Planting
Cost
(euro/acre) |\n", "| :-- | :--: | :--: | :--: | :--: | :--: | :--: |\n", "| Wheat | 170 | - | _ | 238 | 200 | 150 |\n", "| Corn | 150 | - | - | 210 | 240 | 230 |\n", "| Beets | 36 | 6000 | 10 | 0 | - | 260 |\n" ] }, { "cell_type": "markdown", "id": "4b603542-8a50-42dc-8327-412d771d8beb", "metadata": {}, "source": [ "## Data Modeling" ] }, { "cell_type": "code", "execution_count": 3, "id": "d747fa6a-117a-4828-804b-01a0fa7da3ba", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
wheatcornbeets
good3.03.624.0
average2.53.020.0
poor2.02.416.0
\n", "
" ], "text/plain": [ " wheat corn beets\n", "good 3.0 3.6 24.0\n", "average 2.5 3.0 20.0\n", "poor 2.0 2.4 16.0" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "yields = pd.DataFrame(\n", " {\n", " \"good\": {\"wheat\": 3.0, \"corn\": 3.6, \"beets\": 24},\n", " \"average\": {\"wheat\": 2.5, \"corn\": 3.0, \"beets\": 20},\n", " \"poor\": {\"wheat\": 2.0, \"corn\": 2.4, \"beets\": 16},\n", " }\n", ").T\n", "display(yields)" ] }, { "cell_type": "code", "execution_count": 4, "id": "b2e950c8-7a4e-41e1-b0a3-940bed3714c2", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
planting_costsell_pricebuy_costrequireddemandexcess_price
wheat150.0170.0238.0200.020000.00.0
corn230.0150.0210.0240.020000.00.0
beets260.036.020000.00.06000.010.0
\n", "
" ], "text/plain": [ " planting_cost sell_price buy_cost required demand excess_price\n", "wheat 150.0 170.0 238.0 200.0 20000.0 0.0\n", "corn 230.0 150.0 210.0 240.0 20000.0 0.0\n", "beets 260.0 36.0 20000.0 0.0 6000.0 10.0" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "crops = pd.DataFrame(\n", " {\n", " \"wheat\": {\n", " \"planting_cost\": 150,\n", " \"sell_price\": 170,\n", " \"buy_cost\": 238,\n", " \"required\": 200,\n", " },\n", " \"corn\": {\n", " \"planting_cost\": 230,\n", " \"sell_price\": 150,\n", " \"buy_cost\": 210,\n", " \"required\": 240,\n", " },\n", " \"beets\": {\n", " \"planting_cost\": 260,\n", " \"sell_price\": 36,\n", " \"demand\": 6000,\n", " \"excess_price\": 10,\n", " },\n", " }\n", ").T\n", "\n", "# fill NaN\n", "\n", "bigM = 20000\n", "crops[\"buy_cost\"].fillna(bigM, inplace=True)\n", "crops[\"demand\"].fillna(bigM, inplace=True)\n", "crops[\"required\"].fillna(0, inplace=True)\n", "crops[\"excess_price\"].fillna(0, inplace=True)\n", "\n", "display(crops)" ] }, { "cell_type": "markdown", "id": "bdb1ffea-d504-49aa-b985-2974808b0dad", "metadata": {}, "source": [ "## Model Building" ] }, { "cell_type": "code", "execution_count": 5, "id": "631bbb5b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Overwriting farmer.mod\n" ] } ], "source": [ "%%writefile farmer.mod\n", "\n", "set CROPS;\n", "set SCENARIOS;\n", "\n", "param bound;\n", "param total_land;\n", "\n", "param planting_cost{CROPS};\n", "param sell_price{CROPS};\n", "param buy_cost{CROPS};\n", "param required{CROPS};\n", "param demand{CROPS};\n", "param excess_price{CROPS};\n", "\n", "param yields{SCENARIOS, CROPS};\n", "\n", "# first stage variables\n", "var acres{CROPS} >= 0, <= total_land;\n", "\n", "# second stage variables (all units in tons)\n", "var grow{SCENARIOS, CROPS} >= 0, <= bound;\n", "var buy{SCENARIOS, CROPS} >= 0, <= bound;\n", "var sell{SCENARIOS, CROPS} >= 0, <= bound;\n", "var excess{SCENARIOS, CROPS} >= 0, <= bound;\n", "\n", "# first stage model\n", "s.t. limit_on_planted_land: sum{c in CROPS} acres[c] <= total_land;\n", "\n", "# second stage model\n", "s.t. crop_yield {s in SCENARIOS, c in CROPS}:\n", " grow[s, c] == yields[s, c] * acres[c];\n", "s.t. balance {s in SCENARIOS, c in CROPS}:\n", " grow[s, c] + buy[s, c] == sell[s, c] + required[c] + excess[s, c];\n", "\n", "# the selling quota.\n", "# we assume that the excess has smaller sell price,\n", "# so it's optimized out if possible.\n", "s.t. quota {s in SCENARIOS, c in CROPS}:\n", " sell[s, c] <= demand[c];\n", "\n", "# expressions for objective functions\n", "var expense {s in SCENARIOS, c in CROPS} =\n", " planting_cost[c] * acres[c] + buy_cost[c] * buy[s, c];\n", "var revenue {s in SCENARIOS, c in CROPS} =\n", " sell_price[c] * sell[s, c] + excess_price[c] * excess[s, c];\n", "var scenario_profit {s in SCENARIOS} =\n", " sum{c in CROPS}(revenue[s, c] - expense[s, c]);" ] }, { "cell_type": "code", "execution_count": 6, "id": "68d1a2c2-62c8-4d34-a248-9bad1c05ef5b", "metadata": {}, "outputs": [], "source": [ "def farmer(crops, yields, total_land=500):\n", " bound = yields.max().max() * total_land\n", "\n", " m = AMPL()\n", " m.read(\"farmer.mod\")\n", "\n", " # The 'crops' data frame has an AMPL indexing set as index ('CROPS') and AMPL parameters as\n", " # columns ('planting_cost', 'sell_price', 'buy_cost', 'required', 'demand' and 'excess_price').\n", " # It's possible to load all the data in the data frame with the following instruction\n", " m.set_data(crops, \"CROPS\")\n", "\n", " # The 'yields' data frame has an AMPL indexing set as index ('SCENARIOS') and another AMPL\n", " # indexing set as columns ('CROPS').\n", " # We first need to load the data for 'SCENARIOS'.\n", " m.set[\"SCENARIOS\"] = yields.index\n", " # Now we can load the data from the 'yields' data frame.\n", " m.param[\"yields\"] = yields\n", "\n", " m.param[\"bound\"] = bound\n", " m.param[\"total_land\"] = total_land\n", "\n", " m.option[\"solver\"] = SOLVER\n", "\n", " return m" ] }, { "cell_type": "code", "execution_count": 7, "id": "bc71fb08-f582-4836-b8ff-bd4a71289fa0", "metadata": {}, "outputs": [], "source": [ "def farm_report(m):\n", " print(f\"Objective = {m.obj['objective'].value():0.2f}\")\n", "\n", " # indexed by scenarios\n", " df_scenarios = m.var[\"scenario_profit\"].to_pandas()\n", " df_scenarios.index.names = [\"scenario\"]\n", " df_scenarios.rename(columns={df_scenarios.columns[0]: \"profit\"}, inplace=True)\n", "\n", " # indexed by crops\n", " df_crops = m.var[\"acres\"].to_pandas()\n", " df_crops.index.names = [\"crop\"]\n", " df_crops.rename(columns={df_crops.columns[0]: \"acres\"}, inplace=True)\n", "\n", " # indexed by scenarios and crops\n", " df_sc = m.var[\"grow\"].to_pandas()\n", " df_sc.index.names = [\"scenario\", \"crop\"]\n", " df_sc.rename(columns={df_sc.columns[0]: \"grow\"}, inplace=True)\n", "\n", " df_sc[\"buy\"] = m.var[\"buy\"].to_pandas()\n", " df_sc[\"sell\"] = m.var[\"sell\"].to_pandas()\n", " df_sc[\"excess\"] = m.var[\"excess\"].to_pandas()\n", " df_sc[\"expense\"] = m.var[\"expense\"].to_pandas()\n", " df_sc[\"revenue\"] = m.var[\"revenue\"].to_pandas()\n", " df_sc[\"profit\"] = df_sc[\"revenue\"] - df_sc[\"expense\"]\n", "\n", " display(df_scenarios)\n", " display(df_crops)\n", " display(df_sc)" ] }, { "cell_type": "markdown", "id": "6cb1b583-211f-4a2a-beaf-d3a7d869a97b", "metadata": {}, "source": [ "## 1. Mean Solution" ] }, { "cell_type": "code", "execution_count": 8, "id": "ab97bf99-d50c-4c19-9625-853c50b9c6a0", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Objective = 118600.00\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
profit
scenario
mean118600.0
\n", "
" ], "text/plain": [ " profit\n", "scenario \n", "mean 118600.0" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
acres
crop
beets300
corn80
wheat120
\n", "
" ], "text/plain": [ " acres\n", "crop \n", "beets 300\n", "corn 80\n", "wheat 120" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
growbuysellexcessexpenserevenueprofit
scenariocrop
meanbeets6000.006000.0078000216000.0138000.0
corn240.000.00184000.0-18400.0
wheat300.00100.001800017000.0-1000.0
\n", "
" ], "text/plain": [ " grow buy sell excess expense revenue profit\n", "scenario crop \n", "mean beets 6000.0 0 6000.0 0 78000 216000.0 138000.0\n", " corn 240.0 0 0.0 0 18400 0.0 -18400.0\n", " wheat 300.0 0 100.0 0 18000 17000.0 -1000.0" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "m = farmer(crops, pd.DataFrame(yields.mean(), columns=[\"mean\"]).T)\n", "\n", "m.eval(\"maximize objective: sum{s in SCENARIOS} scenario_profit[s];\")\n", "m.get_output(\"solve;\")\n", "\n", "farm_report(m)" ] }, { "cell_type": "markdown", "id": "fd7d91a0-843a-4e9f-8d6c-a38363ab8f32", "metadata": {}, "source": [ "## 2. Stochastic Solution\n" ] }, { "cell_type": "markdown", "id": "86d44583-34a3-49f5-8c91-511c1501731f", "metadata": {}, "source": [ "The problem statement asks for a number of different analyses. In a consulting situation, it is possible the client would ask more \"what if\" questions after hearing the initial results. For these reasons, we build a function that returns an AMPL model and add the variables and expressions needed to address all parts of the problem. " ] }, { "cell_type": "code", "execution_count": 9, "id": "8c8c2d36-4f70-4e1a-a5f2-a45dc0a013ed", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Objective = 108390.00\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
profit
scenario
average109350.0
good167000.0
poor48820.0
\n", "
" ], "text/plain": [ " profit\n", "scenario \n", "average 109350.0\n", "good 167000.0\n", "poor 48820.0" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
acres
crop
beets250
corn80
wheat170
\n", "
" ], "text/plain": [ " acres\n", "crop \n", "beets 250\n", "corn 80\n", "wheat 170" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
growbuysellexcessexpenserevenueprofit
scenariocrop
averagebeets5000.00.05000.0065000.0180000.0115000.0
corn240.00.00.0018400.00.0-18400.0
wheat425.00.0225.0025500.038250.012750.0
goodbeets6000.00.06000.0065000.0216000.0151000.0
corn288.00.048.0018400.07200.0-11200.0
wheat510.00.0310.0025500.052700.027200.0
poorbeets4000.00.04000.0065000.0144000.079000.0
corn192.048.00.0028480.00.0-28480.0
wheat340.00.0140.0025500.023800.0-1700.0
\n", "
" ], "text/plain": [ " grow buy sell excess expense revenue profit\n", "scenario crop \n", "average beets 5000.0 0.0 5000.0 0 65000.0 180000.0 115000.0\n", " corn 240.0 0.0 0.0 0 18400.0 0.0 -18400.0\n", " wheat 425.0 0.0 225.0 0 25500.0 38250.0 12750.0\n", "good beets 6000.0 0.0 6000.0 0 65000.0 216000.0 151000.0\n", " corn 288.0 0.0 48.0 0 18400.0 7200.0 -11200.0\n", " wheat 510.0 0.0 310.0 0 25500.0 52700.0 27200.0\n", "poor beets 4000.0 0.0 4000.0 0 65000.0 144000.0 79000.0\n", " corn 192.0 48.0 0.0 0 28480.0 0.0 -28480.0\n", " wheat 340.0 0.0 140.0 0 25500.0 23800.0 -1700.0" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# maximize mean profit\n", "m = farmer(crops, yields)\n", "\n", "m.eval(\n", " \"maximize objective: (sum{s in SCENARIOS} scenario_profit[s]) / card(SCENARIOS);\"\n", ")\n", "m.get_output(\"solve;\")\n", "\n", "farm_report(m)" ] }, { "cell_type": "markdown", "id": "b724882a-b695-40e9-bcab-753b4110bf48", "metadata": {}, "source": [ "## 3. Worst Case Solution" ] }, { "cell_type": "code", "execution_count": 10, "id": "8eb8181a-fe4b-4e19-adb8-11ee3aad1620", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Objective = 59950.00\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
profit
scenario
average59950.0
good59950.0
poor59950.0
\n", "
" ], "text/plain": [ " profit\n", "scenario \n", "average 59950.0\n", "good 59950.0\n", "poor 59950.0" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
acres
crop
beets375.0
corn25.0
wheat100.0
\n", "
" ], "text/plain": [ " acres\n", "crop \n", "beets 375.0\n", "corn 25.0\n", "wheat 100.0" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
growbuysellexcessexpenserevenueprofit
scenariocrop
averagebeets7500.00.05301.9230772198.07692397500.0212850.0115350.0
corn75.0165.00.0000000.00000040400.00.0-40400.0
wheat250.00.00.00000050.00000015000.00.0-15000.0
goodbeets9000.00.04603.8461544396.15384697500.0209700.0112200.0
corn90.0150.00.0000000.00000037250.00.0-37250.0
wheat300.00.00.000000100.00000015000.00.0-15000.0
poorbeets6000.00.06000.0000000.00000097500.0216000.0118500.0
corn60.0180.00.0000000.00000043550.00.0-43550.0
wheat200.00.00.0000000.00000015000.00.0-15000.0
\n", "
" ], "text/plain": [ " grow buy sell excess expense revenue \n", "scenario crop \n", "average beets 7500.0 0.0 5301.923077 2198.076923 97500.0 212850.0 \\\n", " corn 75.0 165.0 0.000000 0.000000 40400.0 0.0 \n", " wheat 250.0 0.0 0.000000 50.000000 15000.0 0.0 \n", "good beets 9000.0 0.0 4603.846154 4396.153846 97500.0 209700.0 \n", " corn 90.0 150.0 0.000000 0.000000 37250.0 0.0 \n", " wheat 300.0 0.0 0.000000 100.000000 15000.0 0.0 \n", "poor beets 6000.0 0.0 6000.000000 0.000000 97500.0 216000.0 \n", " corn 60.0 180.0 0.000000 0.000000 43550.0 0.0 \n", " wheat 200.0 0.0 0.000000 0.000000 15000.0 0.0 \n", "\n", " profit \n", "scenario crop \n", "average beets 115350.0 \n", " corn -40400.0 \n", " wheat -15000.0 \n", "good beets 112200.0 \n", " corn -37250.0 \n", " wheat -15000.0 \n", "poor beets 118500.0 \n", " corn -43550.0 \n", " wheat -15000.0 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# find worst case profit\n", "m = farmer(crops, yields)\n", "\n", "m.eval(\"var worst_case_profit;\")\n", "m.eval(\n", " \"s.t. lower_bound_profit{s in SCENARIOS}: worst_case_profit <= scenario_profit[s];\"\n", ")\n", "m.eval(\"maximize objective: worst_case_profit;\")\n", "m.get_output(\"solve;\")\n", "\n", "farm_report(m)" ] }, { "cell_type": "code", "execution_count": 11, "id": "73a50577-21cc-476f-a1ad-7b705a2a20f7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Objective = 86600.00\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
profit
scenario
average86600.0
good113250.0
poor59950.0
\n", "
" ], "text/plain": [ " profit\n", "scenario \n", "average 86600.0\n", "good 113250.0\n", "poor 59950.0" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
acres
crop
beets375.0
corn25.0
wheat100.0
\n", "
" ], "text/plain": [ " acres\n", "crop \n", "beets 375.0\n", "corn 25.0\n", "wheat 100.0" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
growbuysellexcessexpenserevenueprofit
scenariocrop
averagebeets7500.00.06000.01500.097500.0231000.0133500.0
corn75.0165.00.00.040400.00.0-40400.0
wheat250.00.050.00.015000.08500.0-6500.0
goodbeets9000.00.06000.03000.097500.0246000.0148500.0
corn90.0150.00.00.037250.00.0-37250.0
wheat300.00.0100.00.015000.017000.02000.0
poorbeets6000.00.06000.00.097500.0216000.0118500.0
corn60.0180.00.00.043550.00.0-43550.0
wheat200.00.00.00.015000.00.0-15000.0
\n", "
" ], "text/plain": [ " grow buy sell excess expense revenue profit\n", "scenario crop \n", "average beets 7500.0 0.0 6000.0 1500.0 97500.0 231000.0 133500.0\n", " corn 75.0 165.0 0.0 0.0 40400.0 0.0 -40400.0\n", " wheat 250.0 0.0 50.0 0.0 15000.0 8500.0 -6500.0\n", "good beets 9000.0 0.0 6000.0 3000.0 97500.0 246000.0 148500.0\n", " corn 90.0 150.0 0.0 0.0 37250.0 0.0 -37250.0\n", " wheat 300.0 0.0 100.0 0.0 15000.0 17000.0 2000.0\n", "poor beets 6000.0 0.0 6000.0 0.0 97500.0 216000.0 118500.0\n", " corn 60.0 180.0 0.0 0.0 43550.0 0.0 -43550.0\n", " wheat 200.0 0.0 0.0 0.0 15000.0 0.0 -15000.0" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "worst_case_profit = m.var[\"worst_case_profit\"].value()\n", "\n", "# maximize expected profit subject to worst case as lower bound\n", "\n", "m = farmer(crops, yields)\n", "\n", "m.eval(\"param worst_case_profit;\")\n", "m.param[\"worst_case_profit\"] = worst_case_profit\n", "m.eval(\n", " \"s.t. lower_bound_profit{s in SCENARIOS}: worst_case_profit <= scenario_profit[s];\"\n", ")\n", "m.eval(\n", " \"maximize objective: (sum{s in SCENARIOS} scenario_profit[s]) / card(SCENARIOS);\"\n", ")\n", "m.get_output(\"solve;\")\n", "\n", "farm_report(m)" ] }, { "cell_type": "markdown", "id": "72d58dac-b8f9-44f9-aebc-1853589f6b98", "metadata": {}, "source": [ "## 4. Risk versus Return" ] }, { "cell_type": "code", "execution_count": 12, "id": "c108ff25-9005-4a3a-99ff-94389fd62816", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAAGwCAYAAACEkkAjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVq0lEQVR4nO3deViU5f4G8Hv2GZZhEdkUEUXFBcylY5hrIqJl0uLeycqy3zlZ2WKdOqaklWVlWXby1DmesrS9zC2UNJcUcUlE3BfMjSVFGJBlBub5/YEzMYLIMDO8A9yf6+LKmfeZmS/fUG/f53mfVyaEECAiIiIip5BLXQARERFRc8JwRUREROREDFdERERETsRwRUREROREDFdERERETsRwRUREROREDFdERERETqSUuoCWxGw248KFC/D29oZMJpO6HCIiIqoHIQSKiooQGhoKufzG56UYrhrRhQsXEBYWJnUZRERE1ABnz55F27ZtbziO4aoReXt7A6j6n6PX6yWuxn2YTCZs2LAB8fHxUKlUUpfTJLGHjmMPHcceOo49dJwremgwGBAWFmb9e/xGGK4akWUqUK/XM1xVYzKZ4OHhAb1ezz9MGog9dBx76Dj20HHsoeNc2cP6LunhgnYiIiIiJ2K4IiIiInIihisiIiIiJ+KaKyIiIjdWWVkJk8kkdRlNhslkglKpRFlZGSorK+v1GpVKBYVC4bQaGK6IiIjckBACOTk5KCgokLqUJkUIgeDgYJw9e9auPSV9fX0RHBzslH0oGa6IiIjckCVYBQYGwsPDg5tP15PZbEZxcTG8vLzqteGnEAIlJSXIy8sDAISEhDhcA8MVERGRm6msrLQGq1atWkldTpNiNpthNBqh1WrrFa4AQKfTAQDy8vIQGBjo8BQhF7QTERG5mYqKCgCAh4eHxJW0HJZeO2N9G8MVERGRmxFCAKj/ppXkOGf2muGKiIiIyIkYroiIiIiciOGKiIiIyIkYrpqBS8XlOHe5BH8UlcNQZkJ5RaV1vp6IiMjdDRkyBDNmzJC6DKfhVgzNwLs/H8dnO3+3eU4mAzRKOTRKRdV/VX/+Wqu6+pzluEpeY6zW+rzt67XVn7vOazVKBRRyLsIkIiLAaDRCrVY32mcpldJHG+krIIdZglR5hdn6nBBAmcmMMpO5jle6jlIuuxq2bhzo1Aog94Ice9cegU6jrD0EXvN6a8ir9j6eGgV0KgWvriGiZkkIgVJT/W7n4mz2/Nk6ZMgQ9OjRA0qlEp9//jmio6Px/vvvY+bMmdi2bRs8PT0RHx+Pd955BwEBAXjggQewZcsWbNmyBYsWLQIAZGVlYfPmzZgxY4bNDvUrV67EXXfdZZ2dSUpKwsqVKzF9+nS8+uqr+P3331FRUQE/Pz/8+9//xk8//YT169ejTZs2ePvtt3HnnXc6vTe1YbhqBuaO6YG5Y3pACAFjpRnlFWaUm8woM1VW/bqi6r9lpkoYK64er6h23Druml9XVKLc9Ofrq/+6tveqMP85FVlhFqgwVuKKsb5/EMixI++Mw72QywBPjRLeGiW8tEp4aZTw0qqqHmuU8NYq4euhgo+HGr46FXx0Kvh6qOCrU8PHo2qcnGfdiMgNlZoq0W32ekk++9DcEfBQ1z8yfPrpp/jb3/6G7du3o6CgALfddhsefvhhvPPOOygtLcXzzz+PcePGYdOmTVi0aBGOHTuGHj16YO7cuQCA1q1b1/uzTpw4ge+++w7ff/+9zeaf8+bNw4IFC/Dmm2/i/fffx+TJk/H777/D39+//t94AzFcNSMymezqWR8FoG38z6+oNFeFO9OfgevPx5W2ge5qcCszVaLUWIEDhw4jPCISJjNsAmG5yYyya0PeNa+3PAcAZgEUlVWgqKwCKLT/e5DLAL1OVRW8rgYwX4+rIay25zxU8NGp4aNTQa3kEkYiIgDo1KkTFixYAAB45ZVX0KtXL7z22mvW40uXLkVYWBiOHTuGzp07Q61Ww8PDA8HBwXZ/ltFoxLJly6yBzGyu+vtgypQpmDhxIgDgtddew3vvvYddu3YhISHB0W/vhhiuyGmUCjmUCjk87JxaN5lMWGc4hFHDO0GlUjXosy2ny4vLKlBUXoHisgoUl1eFrOLyChSXmVBcXgFDWQUKS0woKDWioMSEwtKqr4ISE0pNlTALoKCk6jEuldhVg4daYRPArOHr6pkx25Cmgu/VcR5qTmUS0Y3pVAocmjtCss+2R58+fay/3r9/P3755Rd4eXnVGHfy5El07tzZodrCw8NrPdMVHR1t/bWnpyf0er31/oGuxnBFzYJMJoOHWgkPtRKBDXyP8opKa9AqLDVdDWEmFJQYrc9Xf2x5zlBmghBAibESJcZKXCgss+tzVQoZfKxh7M9gVj2U+XqorGfULGP0OhUvHCBqQSx/zjUFnp6e1l8XFxdj9OjReOONN2qMq+smyXK5vMaV77Xdmqb6Z1V37T/WZTKZ9ayWqzWN/0tEjUCjVCDQW4FAb/vmVCvNAkVlptoDmPWx5SyZ0TqmsMQEY6UZpkqBi8VGXCw2Arhi12frtUr4eKjgo1XBdEWODUUZ8PNSXz1Dpr4a0KoFt6tnz7R2/iuUiKihevfuje+++w7t27e/7pV8arUalZW2a3Rbt26NoqIiXLlyxRqg0tPTXV2uUzBcETlIIZddDS5qhNtx83rLVGb1s2VV/zVWC2l/Pq4+rri86qauhrKqqc6zKAUgx9HCnHp9tlYlv04AqwphPtcs9rdMb3qpueCfiOzz2GOP4eOPP8bEiRPx3HPPwd/fHydOnMCXX36J//znP1AoFGjfvj3S0tJw+vRpeHl5wd/fH/369YOHhwdefPFFPPHEE0hLS8Mnn3wi9bdTLwxXRBKpPpUZ4qOz67WmSrPN1OSlolJs2bkH7Tt3Q3G5+erzRmtAM5T+eUbNbN2moxy5hnK7Plcug/UsmG0Aq1pr5mOduvxz0T8X/BO1bKGhodi+fTuef/55xMfHo7y8HOHh4UhISIBcXvXnwrPPPospU6agW7duKC0tRVZWFtq3b4/PP/8cM2fOxMcff4xhw4YhKSkJ06ZNk/g7ujGGK6ImSKWQI8BLgwAvDQDAZPJC6UmBUbHhdV4UYDYLFBurFvX/OY1pvObs2bWPq/5rWfB/ucSEyyU11z3ciKdaUTOUXb3assZi/2przbh3GVHTsnnz5hrPderUCd9///11X9O5c2ekpqbWeD4xMRGJiYk2zz3yyCPWXyclJSEpKanG6y5fvgy9Xm/zXPX9slyN4YqoBZHLZdBrVdBrVQiz87VlpspqZ8CuE8SuWfB/+YoRReUVEAK4YqzEFWMpzheU2vW5VQv+1X+eIasWwHyuOUvmW+3sGRf8E5FUGK6IqF60KgW0KgUC9Q1b8G9dyH/tFZi1LPa3rDUzVYqrC/7LcbHYvilMoGrBf/WF/NeuI6t+BaaXSoZCI1BuqmzwliBERADDFRG5WPUF//awLPgvKDHh8tUwZrC5IvOaUFZt37JrF/yfya/vpyoxe+9GaFXyWgLYNWfHrFOZfy7499YoOYVJRAxXROSeqi/4D/Vt2IL/ugJY9QX/1sclRgjIUGYyI8dUhhyDfXuWKeR/7llWfbG/r4faJqRdu9bMR6eCSsEF//QnS0i/dp8nch1n9prhioianWsX/NeHyWTCmrXrMOi24SipQLV1ZNfu5v/n9OWfZ9KMKDOZUWkWyL9iRP4Vo901e2mUtSz2V113vRkX/Ddvlv2gSkpKoNPZ948LapiSkqq7cjhjWQDDFRHRVZZ7S7ZSqRBm571dy67uWfbnWjKjTUC79gpMS2grKquawiwur7pVk6ML/qvf7/Laxf7V9zPz1nLBvztTKBTw9fW13q7Fw8ODIbqezGYzjEYjysrKrFs91EUIgZKSEuTl5cHX19fm5s8NxXBFROQElgX/QXYu+K+oNKOorKLWWytVD2I115s5tuBfJgP0WpXt+jEPNXx0Stt7YVZbY2bZKkOj5A7/jcFyE+PGuh9ecyGEQGlpKXQ6nV2B1NfXt0E3jq4NwxURkYSUCjn8PNXw81QDqP0eabURQqDEWFnjTNj11ppZpjEvlxhRYqyEELAGufov+K+iUylqXIHpo1PBW6NA7nkZCnefhb+X1iakccG//WQyGUJCQhAYGFjrPfWodiaTCVu3bsWgQYPqPcWnUqmccsbKguGKiKgJkslk8NQo4amxf8G/scJsDV/Vt8SouuelsdqWGdXXllWNNQug1FSJ0sJKZNd6k3IF1pw5XOvnWhb8/7lZbC1nx66z1kzZghf8KxQKp/7F39wpFApUVFRAq9VKtq0KwxURUQujVsrR2luD1t71X/AP2O7wf7mWKcz84jJkHsuCPiAYhWW248orHFvw761RVl1x6fHnXmU+1Xb29/VQIdhHh4GRAbz/JUmO4YqIiOrFZod/f48ax00mE9aZT2LUqJtqnDEoq3aTcss2GIXXTmdefa76VZqWBf9F5RUoqseC/3F922LBvT2d900TNQDDFRERuVxDF/xXmoV1rZjlTJntVKbReoZs87E/8PWec0joEYzbooJc9J0Q3RjDFRERuS2FXFZtwX/dXl17CB9vy8I/vjuAlKf84ePB2xiRNFruCkEiImpWnonvgg6tPZFXVI6X1xyUuhxqwRiuiIioWdCqFHjz3p6QyYDvfzuPnw/lSl0StVAMV0RE1Gz0CffDwwMiAAAv/nAAhSXcH4oaH8MVERE1K5weJKkxXBERUbNimR6Uc3qQJMJwRUREzU6fcD88PLADAE4PUuNjuCIiombp6eGd/5weXM3pQWo8DFdERNQs2UwP7uP0IDUehisiImq2qk8PvvDDARSU2H9fQyJ7MVwREVGzZpke/KOoHC+vPiR1OdQCMFwREVGzplUp8NbYqunBH/adRwqnB8nFGK6IiKjZ693OD49cnR6c/WMmSowVEldEzRnDFRERtQhPDe+MNr46ZBeWYcnmk1KXQ80YwxUREbUIWpUC/7y9KwDg31tP4Wx+icQVUXPFcEVERC3GyB7BiO3QCuUVZry27rDU5VAzxXBFREQthkwmw5w7u0EuA37KzMGOExelLomaIYYrIiJqUaKC9bjvlnAAwMurD6Gi0ixxRdTcMFwREVGL8/TwzvD1UOFobhGWp52RuhxqZhiuiIioxfH1UOOZ+C4AgIUpx3D5CnduJ+dhuCIiohZp0l/aISrYG4WlJrydclTqcqgZYbgiIqIWSSGXIenO7gCAFWlncOiCQeKKqLlguCIiohbrlg6tcHt0CMwCSFp9EEIIqUuiZoDhioiIWrQXRkVBo5RjV1Y+1h7IlrocagYYroiIqEVr6+eBvw3pCAB4be1hlBorJa6ImjqGKyIiavEeHdQRbXx1uFBYhiVbeN9BcgzDFRERtXg6tQIvjqq67+CSLSdx7jLvO0gNx3BFREQEYFR0MG7p4M/7DpLDJA1XW7duxejRoxEaGgqZTIaVK1faHBdCYPbs2QgJCYFOp0NcXByOHz9uM6Z9+/aQyWQ2X6+//rrNmIyMDAwcOBBarRZhYWFYsGBBjVq++eYbREVFQavVIjo6GuvWrbO7FiIiarpkMhnmjO4OuQxYdyAHO07yvoPUMJKGqytXrqBnz5744IMPaj2+YMECvPfee1iyZAnS0tLg6emJESNGoKyszGbc3LlzkZ2dbf16/PHHrccMBgPi4+MRHh6OvXv34s0330RSUhI++ugj65gdO3Zg4sSJmDp1Kvbt24fExEQkJiYiMzPT7lqIiKjp6hqix+R+VfcdnMv7DlIDKaX88JEjR2LkyJG1HhNC4N1338WsWbMwZswYAMCyZcsQFBSElStXYsKECdax3t7eCA4OrvV9li9fDqPRiKVLl0KtVqN79+5IT0/HwoULMW3aNADAokWLkJCQgJkzZwIA5s2bh5SUFCxevBhLliyxq5bqysvLUV5ebn1sMFRtUGcymWAymexpVbNm6QV70nDsoePYQ8c1lx4+PjQCq/afx5GcInyemoXJ/do12mc3lx5KyRU9tPe9ZMJNdkyTyWT44YcfkJiYCAA4deoUOnbsiH379uGmm26yjhs8eDBuuukmLFq0CEDVtGBZWRlMJhPatWuHSZMm4amnnoJSWZUb77//fhgMBpspx19++QW33XYb8vPz4efnh3bt2uHpp5/GjBkzrGPmzJmDlStXYv/+/fWu5VpJSUl4+eWXazy/YsUKeHh4NKxRRETkcttyZPg2SwEPhcCsXpXwVEldEUmppKQEkyZNQmFhIfR6/Q3HS3rmqi45OTkAgKCgIJvng4KCrMcA4IknnkDv3r3h7++PHTt24IUXXkB2djYWLlxofZ+IiIga72E55ufnh5ycnDo/p761XOuFF17A008/bX1sMBgQFhaG+Pj4ev3PaSlMJhNSUlIwfPhwqFT8E6wh2EPHsYeOa049jK8048CHO3E0txgH5RFIunoloas1px5KxRU9tMw81Zfbhqv6qh5eYmJioFar8eijj2L+/PnQaDQSVgZoNJpaa1CpVPxNUwv2xXHsoePYQ8c1hx6qVMCcO7tj0sdp+GL3WdwX2x5dQxrvH8XNoYdSc2YP7X0ft92KwbKGKjc31+b53Nzc666vAoB+/fqhoqICp0+ftr5Pbe9R/TOuN6b68YbUQkRETVf/jgEYFR0MswBe5n0HyQ5uG64iIiIQHByMjRs3Wp8zGAxIS0tDbGzsdV+Xnp4OuVyOwMBAAEBsbCy2bt1qsxgtJSUFXbp0gZ+fn3VM9c+xjLF8TkNrISKipu3FUV2hUcqx81Q+fsq8/jIQouokDVfFxcVIT09Heno6ACArKwvp6ek4c+YMZDIZZsyYgVdeeQWrVq3CgQMHcP/99yM0NNS66D01NRXvvvuuddH58uXL8dRTT+G+++6zBqdJkyZBrVZj6tSpOHjwIL766issWrTIZjrxySefRHJyMt5++20cOXIESUlJ2LNnD6ZPnw4A9aqFiIian7Z+Hnh0cNV9B1/lfQepniRdc7Vnzx4MHTrU+tgSeKZMmYJPPvkEzz33HK5cuYJp06ahoKAAAwYMQHJyMrRaLYCqNU1ffvklkpKSUF5ejoiICDz11FM2wcnHxwcbNmzAY489hj59+iAgIACzZ8+2bsMAAP3798eKFSswa9YsvPjii+jUqRNWrlyJHj16WMfcqBYiImqe/ja4I77dcxbnC0rx760nMSOus9QlkZuTNFwNGTKkzjlsmUyGuXPnYu7cubUe7927N3bu3HnDz4mJicG2bdvqHDN27FiMHTu2wbUQEVHzpFMr8OLtXTF9xT4s2XISY/uGoY2vTuqyyI257ZorIiIid3F7dAj6RfijzMT7DtKNMVwRERHdQPX7Dq7NyMbOU5ekLoncGMMVERFRPXQL1WPiX6puhZO06iDvO0jXxXBFRERUT8/Ed4Feq8SRnCJ8sfus1OWQm2K4IiIiqid/TzWeie8CAHh7w1EUlvAGy1QTwxUREZEdJvdrh85BXigoMWHxL8elLofcEMMVERGRHZQKOV64eiPnT3f8jrP5JRJXRO6G4YqIiMhOQzq3xoDIABgrzViw/qjU5ZCbYbgiIiKyk0wmwwujoiCTAav3X0D62QKpSyI3wnBFRETUAN1DfXB3r7YAgNfWHq7zjiPUsjBcERERNdCzIzpDo5Rj1+l8bDiUK3U55CYYroiIiBooxEeHhwdGAADe+OkITNxYlMBwRURE5JD/G9wRrTzVOHXxCr7YdUbqcsgNMFwRERE5wFurwozhnQEA7/58HIYybiza0jFcEREROWjCzWHo0NoT+VeMWLL5pNTlkMQYroiIiBykUsjxwsiqjUX/+2sWLhSUSlwRSYnhioiIyAniugbiLxH+KK8w460N3Fi0JWO4IiIicgKZTIZ/Xr0tzg/7ziPzfKHEFZFUGK6IiIicpGeYL8bcFAohgNfWcWPRlorhioiIyImeje8CtUKOHScvYfPRP6QuhyTAcEVEROREYf4eePDW9gCqzl5VcGPRFofhioiIyMn+PjQSvh4qHM8rxo/pF6QuhxoZwxUREZGT+ehUmDaoAwDgP79mce1VC8NwRURE5AKT/tIOOpUCh7MNSD15SepyqBExXBEREbmAr4ca9/ZpC6BqY1FqORiuiIiIXMSysH3jkTyc+qNY2mKo0TBcERERuUiH1l4YFhUIAPjf9tPSFkONhuGKiIjIhaYOjAAAfLv3HApKjBJXQ42B4YqIiMiFYju0QtcQPUpNlVix64zU5VAjYLgiIiJyIZlMhqkDqs5eLdvxO0zcVLTZY7giIiJysdE9Q9DaW4McQxnWHciWuhxyMYYrIiIiF9MoFbj/lnAAwH+2cVPR5o7hioiIqBFMviUcGqUcB84XYvfpy1KXQy7EcEVERNQI/D3VuLu3ZVPRUxJXQ67EcEVERNRIpg5oDwDYcCgXv1+6Im0x5DIMV0RERI0kMtAbgzu3hhDcVLQ5Y7giIiJqRJZtGb7ZcxaGMpPE1ZArMFwRERE1ooGdAtA5yAtXjJX4atdZqcshF2C4IiIiakTVNxX9ZMdpVHBT0WaH4YqIiKiRjbmpDVp5qnG+oBTJB3OkLoecjOGKiIiokWlVCtx3dVPR//6aJXE15GwMV0RERBK475ZwqBVy7DtTgL2/c1PR5oThioiISAKtvTUYc1MoAGApz141KwxXREREEpk6sGph+0+Z2Th3uUTiashZGK6IiIgkEhWsx4DIAJgF8OmO01KXQ07CcEVERCQhy7YMX+46i+LyComrIWdguCIiIpLQ4M6t0aG1J4rKK/Dtb+elLoecgOGKiIhIQnL5n5uKfpp6BmYhcUHkMIYrIiIiid3dqy18PVQ4d7kUB/JlUpdDDmK4IiIikphOrcDkfu0AAJuz+VdzU8f/g0RERG7g/tj2UClkOFUkQ8a5QqnLIQcwXBEREbmBIL0Wt/cIBgD8b8fvEldDjmC4IiIichMP9K+632DywVxkF5ZKXA01FMMVERGRm+geqkekXqDCLPApz141WQxXREREbmRIiBkAsCLtd1zhpqJNEsMVERGRG+nuJxDu7wFDWQW+++2c1OVQAzBcERERuRG5DJgSW7Utw9Jfs2DmrqJNDsMVERGRm7m7Vyj0WiVOXyrBxiN5UpdDdmK4IiIicjOeGiUmXt1U9L+/npK4GrIXwxUREZEbmhLbHgq5DDtP5ePgBW4q2pQwXBEREbmhUF8dRkWHAAD++2uWxNWQPRiuiIiI3NTUAREAgNX7LyDPUCZxNVRfDFdERERu6qYwX/QN94OpUmBZKjcVbSoYroiIiNyY5ezVV3vOoqLSLHE1VB8MV0RERG5sWNcg+Huq8UdROX49cVHqcqgeGK6IiIjcmFopx509QwEA3/92XuJqqD4kDVdbt27F6NGjERoaCplMhpUrV9ocF0Jg9uzZCAkJgU6nQ1xcHI4fP24zJj8/H5MnT4Zer4evry+mTp2K4uJimzEZGRkYOHAgtFotwsLCsGDBghq1fPPNN4iKioJWq0V0dDTWrVtndy1ERESucHfvNgCA9QdzYCgzSVwN3UiDw5XRaMS5c+dw5swZmy97XLlyBT179sQHH3xQ6/EFCxbgvffew5IlS5CWlgZPT0+MGDECZWV/XjExefJkHDx4ECkpKVizZg22bt2KadOmWY8bDAbEx8cjPDwce/fuxZtvvomkpCR89NFH1jE7duzAxIkTMXXqVOzbtw+JiYlITExEZmamXbUQERG5QnQbH0QGeqG8woyfDmRLXQ7dgNLeFxw/fhwPPfQQduzYYfO8EAIymQyVlZX1fq+RI0di5MiRtR4TQuDdd9/FrFmzMGbMGADAsmXLEBQUhJUrV2LChAk4fPgwkpOTsXv3bvTt2xcA8P7772PUqFF46623EBoaiuXLl8NoNGLp0qVQq9Xo3r070tPTsXDhQmsIW7RoERISEjBz5kwAwLx585CSkoLFixdjyZIl9aqlNuXl5SgvL7c+NhgMAACTyQSTif/ysLD0gj1pOPbQceyh49hDx9XVw8SeIXgr5Ti+3XsOd98U0tilNRmu+Dm0973sDlcPPPAAlEol1qxZg5CQEMhkMnvfol6ysrKQk5ODuLg463M+Pj7o168fUlNTMWHCBKSmpsLX19carAAgLi4OcrkcaWlpuOuuu5CamopBgwZBrVZbx4wYMQJvvPEGLl++DD8/P6SmpuLpp5+2+fwRI0ZYpynrU0tt5s+fj5dffrnG8xs2bICHh0eD+tKcpaSkSF1Ck8ceOo49dBx76LjaeuhdDsigwO7Tl/HZ9+vQSitBYU2IM38OS0pK7Bpvd7hKT0/H3r17ERUVZe9L7ZKTkwMACAoKsnk+KCjIeiwnJweBgYE2x5VKJfz9/W3GRERE1HgPyzE/Pz/k5OTc8HNuVEttXnjhBZvQZjAYEBYWhvj4eOj1+jq++5bFZDIhJSUFw4cPh0qlkrqcJok9dBx76Dj20HE36uH6wj3YcTIfhX5d8NehHSWo0P254ufQMvNUX3aHq27duuHiRV4KWh8ajQYajabG8yqVin/w1IJ9cRx76Dj20HHsoeOu18N7+4Rhx8l8/Lg/GzOGd3HZ7FFz4MyfQ3vfx+4F7W+88Qaee+45bN68GZcuXYLBYLD5cpbg4GAAQG5urs3zubm51mPBwcHIy8uzOV5RUYH8/HybMbW9R/XPuN6Y6sdvVAsREZGrjegeDA+1AqcvleC3MwVSl0PXYXe4iouLw86dOzFs2DAEBgbCz88Pfn5+8PX1hZ+fn9MKi4iIQHBwMDZu3Gh9zmAwIC0tDbGxsQCA2NhYFBQUYO/evdYxmzZtgtlsRr9+/axjtm7darMYLSUlBV26dLHWGxsba/M5ljGWz6lPLURERK7mqVEioUfVP+q//+2cxNXQ9dg9LfjLL7847cOLi4tx4sQJ6+OsrCykp6fD398f7dq1w4wZM/DKK6+gU6dOiIiIwEsvvYTQ0FAkJiYCALp27YqEhAQ88sgjWLJkCUwmE6ZPn44JEyYgNLRqw7VJkybh5ZdfxtSpU/H8888jMzMTixYtwjvvvGP93CeffBKDBw/G22+/jdtvvx1ffvkl9uzZY92uQSaT3bAWIiKixnBP77b4/rfzWL3/Al66oxu0KoXUJdG1hB2MRqO47bbbxLFjx+x52XX98ssvAkCNrylTpgghhDCbzeKll14SQUFBQqPRiGHDhomjR4/avMelS5fExIkThZeXl9Dr9eLBBx8URUVFNmP2798vBgwYIDQajWjTpo14/fXXa9Ty9ddfi86dOwu1Wi26d+8u1q5da3O8PrXcSGFhoQAgCgsL7Xpdc2c0GsXKlSuF0WiUupQmiz10HHvoOPbQcfXpYUWlWdzy2s8i/Pk1Ym3GhUasrmlwxc+hvX9/23XmSqVSISMjw2nBbsiQIRBCXPe4TCbD3LlzMXfu3OuO8ff3x4oVK+r8nJiYGGzbtq3OMWPHjsXYsWMdqoWIiMjVFHIZ7urVBv/afBLf/3YOo6K555W7sXvN1X333Yf//ve/rqiFiIiI6sFyO5zNR//AxeLyG4ymxmb3mquKigosXboUP//8M/r06QNPT0+b4wsXLnRacURERFRTZKA3erb1wf5zhVi9/wIevDXixi+iRmN3uMrMzETv3r0BAMeOHbM5xv02iIiIGsfdvdti/7lCfP/beYYrNyPp1YJERETUMKN7hmLemkM4cL4Qx3KL0DnIW+qS6Cq711wRERGR9Pw91RgaVXULuO9/Oy9xNVSd3Weuhg4dWuf036ZNmxwqiIiIiOrnnt5tkHIoFz/sO4eZI7pAIefyHHdgd7i66aabbB6bTCakp6cjMzMTU6ZMcVZdREREdANDowLh66FCrqEcO05exMBOraUuidCAcFV9Z/PqkpKSUFxc7HBBREREVD8apQKjY0Lx2c7f8f1v5xmu3ITT1lzdd999WLp0qbPejoiIiOrBsudVcmYOissrJK6GACeGq9TUVGi1Wme9HREREdXDTWG+6BDgiVJTJZIzc6Quh9CAacG7777b5rEQAtnZ2dizZw9eeuklpxVGRERENyaTyXB37zZ4a8MxfLf3HO7t01bqklo8u89c6fV6+Pj4WL/8/f0xZMgQrFu3DnPmzHFFjURERFSHxF5VU4M7sy4hp7BM4mrI7jNXn3zyiQvKICIiooZq6+eBm9v7Yffpy1h7IBtTB3DHdinZfeaqQ4cOuHTpUo3nCwoK0KFDB6cURURERPa5IyYUALB6/wWJKyG7w9Xp06dRWVlZ4/ny8nKcP88dYomIiKQwMjoYchmQfrYAZ/NLpC6nRav3tOCqVausv16/fj18fHysjysrK7Fx40a0b9/eqcURERFR/QR6axHbsRW2n7iE1RkX8PchkVKX1GLVO1wlJiYCqLoq4dqd2FUqFdq3b4+3337bqcURERFR/Y2OCa0KV/uzGa4kVO9pQbPZDLPZjHbt2iEvL8/62Gw2o7y8HEePHsUdd9zhylqJiIioDgk9gqGUy3A424ATeUVSl9Ni2b3mKisrCwEBAQCAsjJe7klEROQufD3UGNS56hY4q/dnS1xNy2V3uDKbzZg3bx7atGkDLy8vnDp1CgDw0ksv4b///a/TCyQiIqL6G90zBACwOuMChBASV9My2R2uXnnlFXzyySdYsGAB1Gq19fkePXrgP//5j1OLIyIiIvvEdQ2CRinHqT+u4FC2QepyWiS7w9WyZcvw0UcfYfLkyVAoFNbne/bsiSNHjji1OCIiIrKPt1aFoV0CAQBrMjg1KAW7w9X58+cRGVnzCgSz2QyTyeSUooiIiKjhRvf8c0NRTg02PrvDVbdu3bBt27Yaz3/77bfo1auXU4oiIiKihrstKhAeagXOXS5F+tkCqctpcey+t+Ds2bMxZcoUnD9/HmazGd9//z2OHj2KZcuWYc2aNa6okYiIiOygUyswvFsQfky/gNX7s9GrnZ/UJbUodp+5GjNmDFavXo2ff/4Znp6emD17Ng4fPozVq1dj+PDhrqiRiIiI7DT66r0G12RcQKWZU4ONya4zVxUVFXjttdfw0EMPISUlxVU1ERERkYMGdg6AXqtEXlE5dp/Oxy0dWkldUoth15krpVKJBQsWoKKiwlX1EBERkRNolAok9AgGULWwnRqP3dOCw4YNw5YtW1xRCxERETmR5arBnzJzYKo0S1xNy2H3gvaRI0fiH//4Bw4cOIA+ffrA09PT5vidd97ptOKIiIio4WI7tEIrTzUuXTFix8lLGHz11jjkWnaHq7///e8AgIULF9Y4JpPJUFlZ6XhVRERE5DClQo6R0cH4fOcZrN5/geGqkTTo3oLX+2KwIiIici+WqwbXZ+agvIJ/TzcGu8MVERERNR03t/dHkF6DovIKbDn6h9TltAgMV0RERM2YXC7DHVfPXq3mvQYbBcMVERFRM2e5avDnQ7koMXI7JVdjuCIiImrmerb1QZi/DqWmSmw6kid1Oc0ewxUREVEzJ5PJrAvbuaGo69m9FQNQdcXgiRMnkJeXB7PZdlOyQYMGOaUwIiIicp7RPUPxr80n8cuRP2AoM0GvVUldUrNld7jauXMnJk2ahN9//x1C2N4IkvtcERERuaeoYG90bO2Jk39cQcrBXNzTp63UJTVbdk8L/t///R/69u2LzMxM5Ofn4/Lly9av/Px8V9RIREREDpLJZNaF7aszODXoSnafuTp+/Di+/fZbREZGuqIeIiIicpE7YkLx7s/H8evxi8i/YoS/p1rqkpolu89c9evXDydOnHBFLURERORCkYFe6BaiR4VZIDkzR+pymi27z1w9/vjjeOaZZ5CTk4Po6GioVLYL4mJiYpxWHBERETnX6J6hOJRtwOr9FzCpXzupy2mW7A5X99xzDwDgoYcesj4nk8kghOCCdiIiIjd3R0wI3kg+gp1Zl5BnKEOgXit1Sc2O3eEqKyvLFXUQERFRIwjz90Cvdr7Yd6YAaw9k48FbI6QuqdmxO1yFh4e7og4iIiJqJKNjQrHvTAFW77/AcOUCDdpEFAAOHTqEM2fOwGg02jx/5513OlwUERERuc7tMSGYt/YQfjtTgHOXS9DWz0PqkpoVu8PVqVOncNddd+HAgQPWtVZA1borAFxzRURE5OaC9Fr0i/DHzlP5WJuRjUcHd5S6pGbF7q0YnnzySURERCAvLw8eHh44ePAgtm7dir59+2Lz5s0uKJGIiIicjRuKuo7d4So1NRVz585FQEAA5HI55HI5BgwYgPnz5+OJJ55wRY1ERETkZCN7hEAhlyHzvAGn/iiWupxmxe5wVVlZCW9vbwBAQEAALlyoSrzh4eE4evSoc6sjIiIil/D3VGNAZAAAYE1GtsTVNC92h6sePXpg//79AKp2a1+wYAG2b9+OuXPnokOHDk4vkIiIiFzjjpgQAMCq/Resa6jJcXaHq1mzZsFsNgMA5s6di6ysLAwcOBDr1q3De++95/QCiYiIyDXiuwdDrZDjRF4xjuYWSV1Os2H31YIjRoyw/joyMhJHjhxBfn4+/Pz8rFcMEhERkfvz0akwuEtrpBzKxer9FxAVrJe6pGbB7jNXFidOnMD69etRWloKf39/Z9ZEREREjcR61eD+bE4NOond4erSpUsYNmwYOnfujFGjRiE7u2oR3NSpU/HMM884vUAiIiJynbiugdCpFDiTX4KMc4VSl9Ms2B2unnrqKahUKpw5cwYeHn/u6Dp+/HgkJyc7tTgiIiJyLQ+1EsO6BgIAVu/nnlfOYHe42rBhA9544w20bdvW5vlOnTrh999/d1phRERE1DgsU4NrMrJhNnNq0FF2h6srV67YnLGyyM/Ph0ajcUpRRERE1HgGd24Nb40SOYYy7D1zWepymjy7w9XAgQOxbNky62OZTAaz2YwFCxZg6NChTi2OiIiIXE+rUiC+ezAATg06g91bMSxYsADDhg3Dnj17YDQa8dxzz+HgwYPIz8/H9u3bXVEjERERudjoniH47rdzWHcgG7Pv6AalosEbCrR4Ddqh/dixYxgwYADGjBmDK1eu4O6778a+ffvQsSPvqk1ERNQU3RoZAD8PFS4WG7HzVL7U5TRpdp+5AgAfHx/885//dHYtREREJBGVQo6R0SFYkXYGq/dfwIBOAVKX1GQ1KFyVlZUhIyMDeXl51lvhWNx5551OKYyIiIga1+iYUKxIO4OfMrMxL7EH1EpODTaE3eEqOTkZ999/Py5evFjjmEwmQ2VlpVMKIyIiosb1lwh/tPbW4I+icmw99gfiugVJXVKTZHckffzxxzF27FhkZ2fDbDbbfDFYERERNV0KuQy3R4cAANZk8KrBhrI7XOXm5uLpp59GUFDjpNmioiLMmDED4eHh0Ol06N+/P3bv3m09/sADD0Amk9l8JSQk2LxHfn4+Jk+eDL1eD19fX0ydOhXFxcU2YzIyMjBw4EBotVqEhYVhwYIFNWr55ptvEBUVBa1Wi+joaKxbt8413zQREZFELBuKphzKRamRJ00awu5wde+992Lz5s0uKKV2Dz/8MFJSUvDZZ5/hwIEDiI+PR1xcHM6fP28dk5CQgOzsbOvXF198YfMekydPxsGDB5GSkoI1a9Zg69atmDZtmvW4wWBAfHw8wsPDsXfvXrz55ptISkrCRx99ZB2zY8cOTJw4EVOnTsW+ffuQmJiIxMREZGZmur4JREREjaR3O1+08dXhirESvxzNk7qcJsnuNVeLFy/G2LFjsW3bNkRHR0OlUtkcf+KJJ5xWXGlpKb777jv8+OOPGDRoEAAgKSkJq1evxocffohXXnkFAKDRaBAcHFzrexw+fBjJycnYvXs3+vbtCwB4//33MWrUKLz11lsIDQ3F8uXLYTQasXTpUqjVanTv3h3p6elYuHChNYQtWrQICQkJmDlzJgBg3rx5SElJweLFi7FkyZJaP7u8vBzl5eXWxwaDAQBgMplgMpmc0KHmwdIL9qTh2EPHsYeOYw8d5y49HNUjCB//eho/7juH4VFN66pBV/TQ3veyO1x98cUX2LBhA7RaLTZv3gyZTGY9JpPJnBquKioqUFlZCa1Wa/O8TqfDr7/+an28efNmBAYGws/PD7fddhteeeUVtGrVCgCQmpoKX19fa7ACgLi4OMjlcqSlpeGuu+5CamoqBg0aBLVabR0zYsQIvPHGG7h8+TL8/PyQmpqKp59+2qaOESNGYOXKldetf/78+Xj55ZdrPL9hw4ZabyHU0qWkpEhdQpPHHjqOPXQce+g4qXvoewUAlNh0OBffr14HrULSchrEmT0sKSmxa7zd4eqf//wnXn75ZfzjH/+AXO7aSzS9vb0RGxuLefPmoWvXrggKCsIXX3yB1NRUREZGAqiaErz77rsRERGBkydP4sUXX8TIkSORmpoKhUKBnJwcBAYG2ryvUqmEv78/cnJyAAA5OTmIiIiwGWNZU5aTkwM/Pz/k5OTUWGcWFBRkfY/avPDCCzaBzGAwICwsDPHx8dDr9Q1vTDNjMpmQkpKC4cOH1zgTSvXDHjqOPXQce+g4d+mhEALfnt+OrEslULTrhVE9QySrxV6u6KFl5qm+7A5XRqMR48ePd3mwsvjss8/w0EMPoU2bNlAoFOjduzcmTpyIvXv3AgAmTJhgHRsdHY2YmBh07NgRmzdvxrBhwxqlxuvRaDS13sxapVLxD55asC+OYw8dxx46jj10nDv0cPRNbfDexuP4KTMX9/ZtJ2ktDeHMHtr7PnYnpClTpuCrr76y92UN1rFjR2zZsgXFxcU4e/Ysdu3aBZPJhA4dOtQ6vkOHDggICMCJEycAAMHBwcjLs12QV1FRgfz8fOs6reDgYOTm5tqMsTy+0ZjrrfUiIiJqykbHVJ2t2nr8DxSUGCWupmmx+8xVZWUlFixYgPXr1yMmJqZGmlu4cKHTiqvO09MTnp6euHz5MtavX1/rVgkAcO7cOVy6dAkhIVU/FLGxsSgoKMDevXvRp08fAMCmTZtgNpvRr18/65h//vOfMJlM1u8nJSUFXbp0gZ+fn3XMxo0bMWPGDOtnpaSkIDY21iXfLxERkZQ6BXkjKtgbR3KKsP5gDsbf3PTOXknF7jNXBw4cQK9evSCXy5GZmYl9+/ZZv9LT051e4Pr165GcnIysrCykpKRg6NChiIqKwoMPPoji4mLMnDkTO3fuxOnTp7Fx40aMGTMGkZGRGDFiBACga9euSEhIwCOPPIJdu3Zh+/btmD59OiZMmIDQ0Kq9PCZNmgS1Wo2pU6fi4MGD+Oqrr7Bo0SKb9VJPPvkkkpOT8fbbb+PIkSNISkrCnj17MH36dKd/z0RERO7AsufV6v3ZElfStNh95uqXX35xRR3XVVhYiBdeeAHnzp2Dv78/7rnnHrz66qtQqVSoqKhARkYGPv30UxQUFCA0NBTx8fGYN2+ezVqn5cuXY/r06Rg2bBjkcjnuuecevPfee9bjPj4+2LBhAx577DH06dMHAQEBmD17ts1eWP3798eKFSswa9YsvPjii+jUqRNWrlyJHj16NGo/iIiIGsvomFC8uf4odpy8iD+KytHau+Y6YqqpQTdubkzjxo3DuHHjaj2m0+mwfv36G76Hv78/VqxYUeeYmJgYbNu2rc4xY8eOxdixY2/4eURERM1Bu1Ye6Bnmi/1nC/BTZjbuj20vdUlNAm93TURERNdlWdi+ej/vNVhfDFdERER0XbdfDVe7T1/GhYJSiatpGhiuiIiI6LpCfHT4S3t/AMDaDC5srw+GKyIiIqrT6Ks7tK/O4NRgfTBcERERUZ1GRodALgMyzhXi9MUrUpfj9hiuiIiIqE4BXhrcGhkAAFjDs1c3xHBFRERENzQ6pmpD0TVcd3VDDFdERER0QyO6B0OlkOFIThGO5RZJXY5bY7giIiKiG/LxUGFw59YAgDXc86pODFdERERUL9Z7DWZkQwghcTXui+GKiIiI6iWuaxC0KjmyLl7BwQsGqctxWwxXREREVC+eGiWGRQUB4O1w6sJwRURERPV2x9Xb4azJyIbZzKnB2jBcERERUb0NjQqEp1qB8wWl2Hf2stTluCWGKyIiIqo3rUqB+O7BAIBV6ZwarA3DFREREdnFcq/BtQdyUMmpwRoYroiIiMguAyJbw0enwsXicqSduiR1OW6H4YqIiIjsolbKMbJH1dTgat4OpwaGKyIiIrKbZUPRnzKzYao0S1yNe2G4IiIiIrvd0qEVArw0KCgx4dcTF6Uux60wXBEREZHdFHIZbo++OjXIDUVtMFwRERFRg1imBjcczEWZqVLiatwHwxURERE1SO92fgj10aK4vAKbj/4hdTlug+GKiIiIGkQul+H2q7fDWZ3BqUELhisiIiJqMMvU4MbDubhSXiFxNe6B4YqIiIgaLLqND8JbeaDMZMbPh3OlLsctMFwRERFRg8lkMoyOqTp7tYYbigJguCIiIiIHWaYGtxz9A4WlJomrkR7DFRERETmkS7A3Ogd5wVhpxoaDOVKXIzmGKyIiInKYZWqQ9xpkuCIiIiInuOPq1OD2Exdxqbhc4mqkxXBFREREDosI8ER0Gx9UmgV+ymzZU4MMV0REROQUo3te3VC0hd9rkOGKiIiInOL2q+uudp3OR66hTOJqpMNwRURERE7RxleHPuF+EAJY24IXtjNcERERkdOM5r0GGa6IiIjIeUbFhEAuA/adKcDZ/BKpy5EEwxURERE5TaC3Frd0aAWg5d4Oh+GKiIiInMpyO5yWetUgwxURERE5VUL3YCjlMhzKNuBEXrHU5TQ6hisiIiJyKj9PNQZ2CgAArGmBC9sZroiIiMjpLFODq/ZfgBBC4moaF8MVEREROd3wbkFQK+U49ccVHM4ukrqcRsVwRURERE7nrVXhti6BAFrenlcMV0REROQSlqnBNRkta2qQ4YqIiIhc4raoQHioFTibX4r95wqlLqfRMFwRERGRS+jUCsR1DQLQsva8YrgiIiIil6k+NWg2t4ypQYYrIiIicplBnQPgrVUi11CO3afzpS6nUTBcERERkctolAokdA8G0HKuGmS4IiIiIpeyTA2uO5CDikqzxNW4HsMVERERuVT/jq3g76lG/hUjdpy8JHU5LsdwRURERC6lVMgxKvrq1GALuGqQ4YqIiIhcbnRM1dRg8sEclFdUSlyNazFcERERkcvd3N4fQXoNisoqsPXYRanLcSmGKyIiInI5uVyGO66evWruU4MMV0RERNQoLFcN/nw4F6XG5js1yHBFREREjaJnWx+E+etQYqzEpiN5UpfjMgxXRERE1ChkspYxNchwRURERI3GctXgpqN5KCozSVyNazBcERERUaPpGuKNjq09YawwI+VQrtTluATDFRERETUamUxmXdjeXKcGGa6IiIioUVnWXW07fhGXrxglrsb5GK6IiIioUUUGeqFbiB4VZoHkgzlSl+N0DFdERETU6Jrz1CDDFRERETW6O2JCAACppy4hz1AmcTXO5fbhqqioCDNmzEB4eDh0Oh369++P3bt3W48LITB79myEhIRAp9MhLi4Ox48ft3mP/Px8TJ48GXq9Hr6+vpg6dSqKi4ttxmRkZGDgwIHQarUICwvDggULatTyzTffICoqClqtFtHR0Vi3bp1rvmkiIqJmLszfA73a+UIIYN2BbKnLcSq3D1cPP/wwUlJS8Nlnn+HAgQOIj49HXFwczp8/DwBYsGAB3nvvPSxZsgRpaWnw9PTEiBEjUFb2ZwqePHkyDh48iJSUFKxZswZbt27FtGnTrMcNBgPi4+MRHh6OvXv34s0330RSUhI++ugj65gdO3Zg4sSJmDp1Kvbt24fExEQkJiYiMzOz8ZpBRETUjFj2vFrVzKYGlVIXUJfS0lJ89913+PHHHzFo0CAAQFJSElavXo0PP/wQ8+bNw7vvvotZs2ZhzJgxAIBly5YhKCgIK1euxIQJE3D48GEkJydj9+7d6Nu3LwDg/fffx6hRo/DWW28hNDQUy5cvh9FoxNKlS6FWq9G9e3ekp6dj4cKF1hC2aNEiJCQkYObMmQCAefPmISUlBYsXL8aSJUtqrb+8vBzl5eXWxwaDAQBgMplgMjXPjdMawtIL9qTh2EPHsYeOYw8d19J6GN81APPWAr+dKcDvfxgQ6qtz+D1d0UN738utw1VFRQUqKyuh1WptntfpdPj111+RlZWFnJwcxMXFWY/5+PigX79+SE1NxYQJE5CamgpfX19rsAKAuLg4yOVypKWl4a677kJqaioGDRoEtVptHTNixAi88cYbuHz5Mvz8/JCamoqnn37apo4RI0Zg5cqV161//vz5ePnll2s8v2HDBnh4eNjbjmYvJSVF6hKaPPbQceyh49hDx7WkHnb0luOEQY6F327GbaHCae/rzB6WlJTYNd6tw5W3tzdiY2Mxb948dO3aFUFBQfjiiy+QmpqKyMhI5ORUXb4ZFBRk87qgoCDrsZycHAQGBtocVyqV8Pf3txkTERFR4z0sx/z8/JCTk1Pn59TmhRdesAlkBoMBYWFhiI+Ph16vt6cVzZrJZEJKSgqGDx8OlUoldTlNEnvoOPbQceyh41piDwsCzmLO6sM4afLDW6Nucfj9XNFDy8xTfbl1uAKAzz77DA899BDatGkDhUKB3r17Y+LEidi7d6/Upd2QRqOBRqOp8bxKpWoxv2nswb44jj10HHvoOPbQcS2ph3f0bIO5a48g84IB5wqNiAjwdMr7OrOH9r6P2y9o79ixI7Zs2YLi4mKcPXsWu3btgslkQocOHRAcHAwAyM21vTdRbm6u9VhwcDDy8vJsjldUVCA/P99mTG3vYTlW1xjLcSIiIrJfKy8Nbo0MAACsaSYL290+XFl4enoiJCQEly9fxvr16zFmzBhEREQgODgYGzdutI4zGAxIS0tDbGwsACA2NhYFBQU2Z7o2bdoEs9mMfv36Wcds3brVZsFaSkoKunTpAj8/P+uY6p9jGWP5HCIiImqY0Vf3vFqdwXDVKNavX4/k5GRkZWUhJSUFQ4cORVRUFB588EHIZDLMmDEDr7zyClatWoUDBw7g/vvvR2hoKBITEwEAXbt2RUJCAh555BHs2rUL27dvx/Tp0zFhwgSEhlZdAjpp0iSo1WpMnToVBw8exFdffYVFixbZrJd68sknkZycjLfffhtHjhxBUlIS9uzZg+nTp0vRFiIiomYjvnsw1Ao5juUW42hOkdTlOMztw1VhYSEee+wxREVF4f7778eAAQOwfv166/znc889h8cffxzTpk3DzTffjOLiYiQnJ9tcYbh8+XJERUVh2LBhGDVqFAYMGGCzh5WPjw82bNiArKws9OnTB8888wxmz55tsxdW//79sWLFCnz00Ufo2bMnvv32W6xcuRI9evRovGYQERE1Qz46FQZ3aQ2gedwOx+0XtI8bNw7jxo277nGZTIa5c+di7ty51x3j7++PFStW1Pk5MTEx2LZtW51jxo4di7Fjx9ZdMBEREdltdM9QpBzKxeqMC3gmvjNkMpnUJTWY25+5IiIiouYvrmsgdCoFfr9UggPnC6UuxyEMV0RERCQ5D7USw7pW7UvZ1KcGGa6IiIjILdxx9V6DazKyYTY7b7f2xsZwRURERG5hSJfW8NIokV1Yhr1nLktdToMxXBEREZFb0KoUiO9edau5pryhKMMVERERuY3RPaumBtceyEZFpVniahqG4YqIiIjcxoDIAPh6qHCx2Ii0rHypy2kQhisiIiJyGyqFHCN7XL0dThOdGmS4IiIiIrcyumdVuPopMwfGiqY3NchwRURERG6lX0QrtPbWoLDUhF9P/CF1OXZjuCIiIiK3opDLcHu0ZWowW+Jq7MdwRURERG7HMjW44WAOykyVEldjH4YrIiIicju9wvzQxleHK8ZK/HIkT+py7MJwRURERG5HLpfhjpirU4MZTeuqQYYrIiIickuWDUU3Hs5DcXmFxNXUH8MVERERuaXuoXpEBHiivMKMnw/lSl1OvTFcERERkVuSyWQYHdP0NhRluCIiIiK3ZZka3Hr8DxSWmCSupn4YroiIiMhtdQryRlSwN0yVAusP5khdTr0wXBEREZFba2pXDTJcERERkVu7I6ZqanD7iYu4WFwucTU3xnBFREREbq19gCdi2vrALICfDrj/7XAYroiIiMjtjb569qop3GuQ4YqIiIjc3u1X113tOp2P7MJSiaupG8MVERERub1QXx1ubu8HAFib4d5nrxiuiIiIqEmw7Hnl7huKMlwRERFRkzCyRwjkMmD/uUKcuVQidTnXxXBFRERETUJrbw36dwwA4N57XjFcERERUZMxuqf732uQ4YqIiIiajBHdg6FSyHAkpwjHc4ukLqdWDFdERETUZPh6qDGoU2sAwGo3vWqQ4YqIiIiaFMtVg2v2X4AQQuJqamK4IiIioiYlrlsQNEo5Tl28goMXDFKXUwPDFRERETUpXholhnUNBOCeVw0yXBEREVGTY7nX4Jr92W43NchwRURERE3O0KhAeKoVOF9Qin1nC6QuxwbDFRERETU5WpUCw7sFAXC/Pa8YroiIiKhJslw1uDYjG5Vm95kaZLgiIiKiJmlgp9bQa5XIKyrHrqx8qcuxYrgiIiKiJkmtlGNkj6u3w3GjqwYZroiIiKjJskwN/nQgG6ZKs8TVVGG4IiIioibrlg7+CPBS43KJCdtPXJS6HAAMV0RERNSEKRVyjIq+OjW43z3uNchwRURERE2aZWpww8EclJsqJa6G4YqIiIiauD7t/BDio0VReQW2nbgkdTkMV0RERNS0yeUy3BETAi+NErlF5VKXw3BFRERETd9jQyOxZ1YcJv8lTOpSoJS6ACIiIiJH+XqoAQAmk/TbMfDMFREREZETMVwRERERORHDFREREZETMVwRERERORHDFREREZETMVwRERERORHDFREREZETMVwRERERORHDFREREZETMVwRERERORHDFREREZETMVwRERERORHDFREREZETKaUuoCURQgAADAaDxJW4F5PJhJKSEhgMBqhUKqnLaZLYQ8exh45jDx3HHjrOFT20/L1t+Xv8RhiuGlFRUREAICwsTOJKiIiIyF5FRUXw8fG54TiZqG8MI4eZzWZcuHAB3t7ekMlkUpfjNgwGA8LCwnD27Fno9Xqpy2mS2EPHsYeOYw8dxx46zhU9FEKgqKgIoaGhkMtvvKKKZ64akVwuR9u2baUuw23p9Xr+YeIg9tBx7KHj2EPHsYeOc3YP63PGyoIL2omIiIiciOGKiIiIyIkYrkhyGo0Gc+bMgUajkbqUJos9dBx76Dj20HHsoePcoYdc0E5ERETkRDxzRUREROREDFdERERETsRwRUREROREDFdERERETsRwRXZLSkqCTCaz+YqKirIeLysrw2OPPYZWrVrBy8sL99xzD3Jzc23e48yZM7j99tvh4eGBwMBAzJw5ExUVFTZjNm/ejN69e0Oj0SAyMhKffPJJjVo++OADtG/fHlqtFv369cOuXbtc8j07W109zM/Px+OPP44uXbpAp9OhXbt2eOKJJ1BYWGjzHuxh3T+HFkIIjBw5EjKZDCtXrrQ5xh7euIepqam47bbb4OnpCb1ej0GDBqG0tNR6PD8/H5MnT4Zer4evry+mTp2K4uJim/fIyMjAwIEDodVqERYWhgULFtSo5ZtvvkFUVBS0Wi2io6Oxbt0613zTTnajHubk5OCvf/0rgoOD4enpid69e+O7776zeY+W3kMAOH/+PO677z60atUKOp0O0dHR2LNnj/W4EAKzZ89GSEgIdDod4uLicPz4cZv3cKs+CiI7zZkzR3Tv3l1kZ2dbv/744w/r8f/7v/8TYWFhYuPGjWLPnj3illtuEf3797cer6ioED169BBxcXFi3759Yt26dSIgIEC88MIL1jGnTp0SHh4e4umnnxaHDh0S77//vlAoFCI5Odk65ssvvxRqtVosXbpUHDx4UDzyyCPC19dX5ObmNk4jHFBXDw8cOCDuvvtusWrVKnHixAmxceNG0alTJ3HPPfdYX88e3vjn0GLhwoVi5MiRAoD44YcfrM+zhzfu4Y4dO4Rerxfz588XmZmZ4siRI+Krr74SZWVl1jEJCQmiZ8+eYufOnWLbtm0iMjJSTJw40Xq8sLBQBAUFicmTJ4vMzEzxxRdfCJ1OJ/79739bx2zfvl0oFAqxYMECcejQITFr1iyhUqnEgQMHGqcRDrhRD4cPHy5uvvlmkZaWJk6ePCnmzZsn5HK5+O2336xjWnoP8/PzRXh4uHjggQdEWlqaOHXqlFi/fr04ceKEdczrr78ufHx8xMqVK8X+/fvFnXfeKSIiIkRpaal1jDv1keGK7DZnzhzRs2fPWo8VFBQIlUolvvnmG+tzhw8fFgBEamqqEEKIdevWCblcLnJycqxjPvzwQ6HX60V5ebkQQojnnntOdO/e3ea9x48fL0aMGGF9/Je//EU89thj1seVlZUiNDRUzJ8/3+Hv0dXq6mFtvv76a6FWq4XJZBJCsIdC1K+H+/btE23atBHZ2dk1whV7eOMe9uvXT8yaNeu6xw8dOiQAiN27d1uf++mnn4RMJhPnz58XQgjxr3/9S/j5+Vl7KoQQzz//vOjSpYv18bhx48Ttt99e47MfffRRe7+lRnejHnp6eoply5bZPOfv7y8+/vhjIQR7KETV9zJgwIDrHjebzSI4OFi8+eab1ucKCgqERqMRX3zxhRDC/frIaUFqkOPHjyM0NBQdOnTA5MmTcebMGQDA3r17YTKZEBcXZx0bFRWFdu3aITU1FUDVNEN0dDSCgoKsY0aMGAGDwYCDBw9ax1R/D8sYy3sYjUbs3bvXZoxcLkdcXJx1jLu7Xg9rU1hYCL1eD6Wy6nag7GGVunpYUlKCSZMm4YMPPkBwcHCN17KHVa7Xw7y8PKSlpSEwMBD9+/dHUFAQBg8ejF9//dX62tTUVPj6+qJv377W5+Li4iCXy5GWlmYdM2jQIKjVauuYESNG4OjRo7h8+bJ1TF19dnd1/Rz2798fX331FfLz82E2m/Hll1+irKwMQ4YMAcAeAsCqVavQt29fjB07FoGBgejVqxc+/vhj6/GsrCzk5OTYfH8+Pj7o16+fzd8r7tRHhiuyW79+/fDJJ58gOTkZH374IbKysjBw4EAUFRUhJycHarUavr6+Nq8JCgpCTk4OgKo1CNX/QrMctxyra4zBYEBpaSkuXryIysrKWsdY3sOd1dXDa128eBHz5s3DtGnTrM+xhzfu4VNPPYX+/ftjzJgxtb6ePay7h6dOnQJQtabokUceQXJyMnr37o1hw4ZZ17rk5OQgMDDQ5j2VSiX8/f2d8vu9qfcQAL7++muYTCa0atUKGo0Gjz76KH744QdERkYCYA8B4NSpU/jwww/RqVMnrF+/Hn/729/wxBNP4NNPPwXw5/dY1/fnbn1U2jWaCMDIkSOtv46JiUG/fv0QHh6Or7/+GjqdTsLKmo66ejh16lTrMYPBgNtvvx3dunVDUlKSBJW6r7p62Lp1a2zatAn79u2TsEL3V1cPu3btCgB49NFH8eCDDwIAevXqhY0bN2Lp0qWYP3++JDW7mxv9Xn7ppZdQUFCAn3/+GQEBAVi5ciXGjRuHbdu2ITo6WsLK3YfZbEbfvn3x2muvAaj6OcvMzMSSJUswZcoUiatrGJ65Iof5+vqic+fOOHHiBIKDg2E0GlFQUGAzJjc31zo1ExwcXOPqQcvjG43R6/XQ6XQICAiAQqGodUxtU0DurnoPLYqKipCQkABvb2/88MMPUKlU1mPsYU3Ve7hp0yacPHkSvr6+UCqV1unUe+65xzodwx7WVL2HISEhAIBu3brZjOnatat12is4OBh5eXk2xysqKpCfn++U3+9NvYcnT57E4sWLsXTpUgwbNgw9e/bEnDlz0LdvX3zwwQcA2EMACAkJueHPGYA6vz936yPDFTmsuLgYJ0+eREhICPr06QOVSoWNGzdajx89ehRnzpxBbGwsACA2NhYHDhyw+Y2QkpICvV5v/Q0WGxtr8x6WMZb3UKvV6NOnj80Ys9mMjRs3Wsc0JdV7CFSdsYqPj4darcaqVaug1WptxrOHNVXv4T/+8Q9kZGQgPT3d+gUA77zzDv73v/8BYA9rU72H7du3R2hoKI4ePWoz5tixYwgPDwdQ1Z+CggLs3bvXenzTpk0wm83o16+fdczWrVthMpmsY1JSUtClSxf4+flZx9TV56akeg9LSkoAVK3Dq06hUMBsNgNgDwHg1ltvrfPnLCIiAsHBwTbfn8FgQFpams3fK27VR7uWvxMJIZ555hmxefNmkZWVJbZv3y7i4uJEQECAyMvLE0JUbcXQrl07sWnTJrFnzx4RGxsrYmNjra+3XAIfHx8v0tPTRXJysmjdunWtl8DPnDlTHD58WHzwwQe1XgKv0WjEJ598Ig4dOiSmTZsmfH19ba7+cld19bCwsFD069dPREdHixMnTthc4l1RUSGEYA+FuPHP4bVwna0Y2MPr9/Cdd94Rer1efPPNN+L48eNi1qxZQqvV2lwin5CQIHr16iXS0tLEr7/+Kjp16mRz+XtBQYEICgoSf/3rX0VmZqb48ssvhYeHR43L35VKpXjrrbfE4cOHxZw5c5rMNgJ19dBoNIrIyEgxcOBAkZaWJk6cOCHeeustIZPJxNq1a63v0dJ7uGvXLqFUKsWrr74qjh8/LpYvXy48PDzE559/bh3z+uuvC19fX/Hjjz+KjIwMMWbMmFq3YnCXPjJckd3Gjx8vQkJChFqtFm3atBHjx4+3+cO2tLRU/P3vfxd+fn7Cw8ND3HXXXSI7O9vmPU6fPi1GjhwpdDqdCAgIEM8884x1mwGLX375Rdx0001CrVaLDh06iP/97381ann//fdFu3bthFqtFn/5y1/Ezp07XfI9O1tdPfzll18EgFq/srKyrO/BHtb9c3ita8OVEOxhfXo4f/580bZtW+Hh4SFiY2PFtm3bbI5funRJTJw4UXh5eQm9Xi8efPBBUVRUZDNm//79YsCAAUKj0Yg2bdqI119/vUYtX3/9tejcubNQq9Wie/fuNuHDnd2oh8eOHRN33323CAwMFB4eHiImJqbG1gwtvYdCCLF69WrRo0cPodFoRFRUlPjoo49sjpvNZvHSSy+JoKAgodFoxLBhw8TRo0dtxrhTH2VCCGHfuS4iIiIiuh6uuSIiIiJyIoYrIiIiIidiuCIiIiJyIoYrIiIiIidiuCIiIiJyIoYrIiIiIidiuCIiIiJyIoYrIiIiIidiuCIiIiJyIoYrIiIHyGQyrFy5UuoyiMiNMFwREV2H0WiUugQiaoIYroioSVqzZg18fX1RWVkJAEhPT4dMJsM//vEP65iHH34Y9913n/Xxd999h+7du0Oj0aB9+/Z4++23bd6zffv2mDdvHu6//37o9XpMmzYNRqMR06dPR0hICLRaLcLDwzF//nzreAC46667IJPJrI9rc+7cOUycOBH+/v7w9PRE3759kZaWBgA4efIkxowZg6CgIHh5eeHmm2/Gzz//bPP6f/3rX+jUqRO0Wi2CgoJw7733Wo+ZzWbMnz8fERER0Ol06NmzJ7799lv7m0pETqGUugAiooYYOHAgioqKsG/fPvTt2xdbtmxBQEAANm/ebB2zZcsWPP/88wCAvXv3Yty4cUhKSsL48eOxY8cO/P3vf0erVq3wwAMPWF/z1ltvYfbs2ZgzZw4A4L333sOqVavw9ddfo127djh79izOnj0LANi9ezcCAwPxv//9DwkJCVAoFLXWWlxcjMGDB6NNmzZYtWoVgoOD8dtvv8FsNluPjxo1Cq+++io0Gg2WLVuG0aNH4+jRo2jXrh327NmDJ554Ap999hn69++P/Px8bNu2zfr+8+fPx+eff44lS5agU6dO2Lp1K+677z60bt0agwcPdmbbiag+BBFRE9W7d2/x5ptvCiGESExMFK+++qpQq9WiqKhInDt3TgAQx44dE0IIMWnSJDF8+HCb18+cOVN069bN+jg8PFwkJibajHn88cfFbbfdJsxmc601ABA//PBDnXX++9//Ft7e3uLSpUv1/t66d+8u3n//fSGEEN99953Q6/XCYDDUGFdWViY8PDzEjh07bJ6fOnWqmDhxYr0/j4ich9OCRNRkDR48GJs3b4YQAtu2bcPdd9+Nrl274tdff8WWLVsQGhqKTp06AQAOHz6MW2+91eb1t956K44fP26dWgSAvn372ox54IEHkJ6eji5duuCJJ57Ahg0b7K4zPT0dvXr1gr+/f63Hi4uL8eyzz6Jr167w9fWFl5cXDh8+jDNnzgAAhg8fjvDwcHTo0AF//etfsXz5cpSUlAAATpw4gZKSEgwfPhxeXl7Wr2XLluHkyZN210pEjuO0IBE1WUOGDMHSpUuxf/9+qFQqREVFYciQIdi8eTMuX77coCkxT09Pm8e9e/dGVlYWfvrpJ/z8888YN24c4uLi7FrTpNPp6jz+7LPPIiUlBW+99RYiIyOh0+lw7733WhfUe3t747fffsPmzZuxYcMGzJ49G0lJSdi9ezeKi4sBAGvXrkWbNm1s3lej0dS7RiJyHoYrImqyLOuu3nnnHWuQGjJkCF5//XVcvnwZzzzzjHVs165dsX37dpvXb9++HZ07d77uWikLvV6P8ePHY/z48bj33nuRkJCA/Px8+Pv7Q6VS2Zz5qk1MTAz+85//WF9zre3bt+OBBx7AXXfdBaDqTNbp06dtxiiVSsTFxSEuLg5z5syBr68vNm3ahOHDh0Oj0eDMmTNcX0XkJhiuiKjJ8vPzQ0xMDJYvX47FixcDAAYNGoRx48bBZDLZhI1nnnkGN998M+bNm4fx48cjNTUVixcvxr/+9a86P2PhwoUICQlBr169IJfL8c033yA4OBi+vr4Aqq4Y3LhxI2699VZoNBr4+fnVeI+JEyfitddeQ2JiIubPn4+QkBDs27cPoaGhiI2NRadOnfD9999j9OjRkMlkeOmll6yL3YGqKyNPnTqFQYMGwc/PD+vWrYPZbEaXLl3g7e2NZ599Fk899RTMZjMGDBiAwsJCbN++HXq9HlOmTHFCp4nILlIv+iIicsSTTz4pAIjDhw9bn+vZs6cIDg6uMfbbb78V3bp1EyqVSrRr1866GN4iPDxcvPPOOzbPffTRR+Kmm24Snp6eQq/Xi2HDhonffvvNenzVqlUiMjJSKJVKER4eft06T58+Le655x6h1+uFh4eH6Nu3r0hLSxNCCJGVlSWGDh0qdDqdCAsLE4sXLxaDBw8WTz75pBBCiG3btonBgwcLPz8/odPpRExMjPjqq6+s7202m8W7774runTpIlQqlWjdurUYMWKE2LJlS33bSEROJBNCCKkDHhEREVFzwasFiYiIiJyI4YqIiIjIiRiuiIiIiJyI4YqIiIjIiRiuiIiIiJyI4YqIiIjIiRiuiIiIiJyI4YqIiIjIiRiuiIiIiJyI4YqIiIjIiRiuiIiIiJzo/wFPGRV6cBqYOAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df = pd.DataFrame()\n", "\n", "for wc in np.linspace(48820, 59950, 50):\n", " m = farmer(crops, yields)\n", "\n", " m.eval(\"param wc;\")\n", " m.param[\"wc\"] = wc\n", " m.eval(\"s.t. min_profit{s in SCENARIOS}: scenario_profit[s] >= wc;\")\n", " # maximize mean profit\n", " m.eval(\n", " \"maximize objective: (sum{s in SCENARIOS} scenario_profit[s]) / card(SCENARIOS);\"\n", " )\n", " m.get_output(\"solve;\")\n", "\n", " df.loc[wc, \"return\"] = m.obj[\"objective\"].value()\n", " df.loc[wc, \"wc\"] = wc\n", "\n", "df.plot(x=\"wc\", y=\"return\", xlabel=\"worst case\", ylabel=\"mean return\", grid=True)" ] }, { "cell_type": "markdown", "id": "f40722ea-5cd9-4207-9877-09b5bf040098", "metadata": {}, "source": [ "## 5. Recommended Planting Strategy" ] }, { "cell_type": "code", "execution_count": 13, "id": "fe026b43-3e02-4de4-9c02-b677fc61a3e5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAGwCAYAAABrUCsdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABdBElEQVR4nO3deVxU9f4/8NeZnW3YdxEwd8slLIXUNFGym2bLL5cS6mbe9szK6n7NrcUblq2Wt3ItLS3LLL0qaqQp4oqZmisiKpsiDvsMzOf3BzI5gsjADAeY1/PxmIfOOWc+5/05HvTlOZ/5HEkIIUBEREREjaaQuwAiIiKi1oLBioiIiMhOGKyIiIiI7ITBioiIiMhOGKyIiIiI7ITBioiIiMhOGKyIiIiI7EQldwHOxGw249y5c/Dw8IAkSXKXQ0RERPUghEBhYSFCQkKgUNR9TYrBqgmdO3cOYWFhcpdBREREDZCZmYk2bdrUuQ2DVRPy8PAAUPUHo9fr7dq2yWTChg0bMHToUKjVaru23RKw/87df4DHwNn7D/AYsP+O67/BYEBYWJjl3/G6MFg1oerbf3q93iHBytXVFXq93ml/oNh/5+0/wGPg7P0HeAzYf8f3vz7DeDh4nYiIiMhOGKyIiIiI7ITBioiIiMhOOMaKiIioGTObzTAajdfdzmQyQaVSoaysDJWVlU1QWfPSmP6r1WoolUq71MFgRURE1EwZjUakp6fDbDZfd1shBIKCgpCZmemUcyU2tv9eXl4ICgpq9LFjsCIiImqGhBDIysqCUqlEWFjYdSemNJvNKCoqgru7+3W3bY0a2n8hBEpKSpCbmwsACA4OblQdDFZERETNUEVFBUpKShASEgJXV9frbl99y1Cn0zltsGpo/11cXAAAubm5CAgIaNRtQec78kRERC1A9TghjUYjcyXOoTq8mkymRrXDYEVERNSMOeN4KTnY6zgzWBERERHZCYMVERERkZ0wWBERERHZCYNVK5BfbERGfgmMzjcfHBERtQIDBw7ExIkT5S7DLjjdQiuwat9ZzPzlEAAVZuzfjAAPLfw9tAjQ6+DvrkWAXmv5NcBDB38PLbxd1RwQSUREDmc0Gpvsm431maHe0RisWoHyCjN0agXKTGYUlVegqLwCJ88X1/kZtVKCn/vlAHY5iPlfDl2WYOahhZ+7Fjq1fab5JyKihhNCoNR07VsTZrMZpcZKqIwVdp/HykWtrPd/xgcOHIgbb7wRKpUKX3/9NW666SZ8/PHHePnll7F161a4ublh6NCheP/99+Hn54dHHnkEv/32G3777Td8+OGHAID09HQkJydj4sSJKCgosLS9atUq3HvvvRBCAACmT5+OVatW4ZlnnsFbb72FjIwMVFRUQJIkfPHFF1izZg3Wr1+P0NBQvPfeexgxYoRdj0ttGKxagScH3oDHYsLw48//Q8/o23GxrBK5heXIKyxHbmEZ8qp/byhHXlE58ouNMFUKZF0qQ9alsuu27+mitgStv0PX1SFMB72LilfBiIgcpNRUia5T18uy70Mz4+CqqX9kWLx4MZ588kls27YNBQUFuOOOOzB+/Hi8//77KC0txSuvvIIHH3wQmzdvxocffoijR4/ixhtvxMyZMwEA/v7+9d7X8ePHsXLlSnz//fcoLS21LJ8xYwYSExMxe/ZsfPzxx3jooYeQkZEBHx+f+ne8ARisWglJkqBTAe383aBWq+vc1lhhxoXiy0GrsNwqhFX/vvplrDTjUqkJl0pNOJ5bVGe7GpWixq1Hf3ddjVuRvu4aqJUc3kdE1Fp16NABiYmJAIA333wTvXr1wttvv21Zv2DBAoSFheHo0aPo2LEjNBoNXF1dERQUZPO+jEYjlixZAl9fXxgMBsvyRx55BGPGjAEAvP322/joo4+wc+dO3HnnnY3sXd0YrJyQRqVAsKcLgj1d6txOCIFLpaYa4av6ffUVsFxDGQxlFTBWmHG2oBRnC0rrbFeSAB9XzeXbjzWvgF05RsxNU//Lz0RErZmLWolDM+Ouud5sNqPQUAgPvYdDbgXaIioqyvL7/fv349dff4W7u3uN7U6cOIGOHTs2qrbw8HD4+/vXeFB19+7dLb93c3ODXq+3PA/QkRis6JokSYKXqwZerhp0CPSoc9syU2XVVa6i6ith1le/LL8vKkelWeBCsREXio34K7uwznZd1MqrroBdHpTvYT0+TK/hFTAiat0kSarzdpzZbEaFRglXjUr2ZwW6ublZfl9UVIThw4fjnXfeqbFdXQ88VigUlrFU1Wp73MyV+7rS1XdvJEmqEb4cgcGK7EKnViLMxxVhPnU/KNRsFsgvMV5x1asMeUVXha/Ly4uNlSg1VSLjQgkyLpTU2a5CAtxVSvz3VAoC9LrLV75qjgPz99DCRcPB+ERETeXmm2/GypUrERERAZWq9tih0Wgsz0as5u/vj8LCQhQXF1vCU1pamqPLbTQGK2pSCkXVtxH93LXocu3/qAAAissrrgpcZVfckvx7+YXicpgFYDBJOJRViENZdV8F89Cq4H/l1a+rx4VdDmGckoKIqPGefvppfPHFFxgzZgwmT54MHx8fHD9+HN9++y2+/PJLKJVKREREIDU1FadOnYK7uzt8fHzQp08fuLq64t///jeee+45pKamYtGiRXJ357oYrKjZctOq4KZVIcKv9su81SoqzcguKMaqdZvRqectyC+pqDkOrLAMuYZylFeYUVhegcK8CpzMq9+UFHVNRxGg18HPXQOtilfBiIhqExISgm3btuGVV17B0KFDUV5ejvDwcNx5552WW5YvvfQSEhIS0LVrV5SWliI9PR0RERH4+uuv8fLLL+OLL77A4MGDMX36dEyYMEHmHtWNwYpaPJVSgUC9DmHuwMCO/tf8VqQQAoXVV8GuGHh/9Tiw3MIyXCwx2TQlhZeruuY4MKsrYFXBTK/jlBRE1LolJyfXWNahQwf88MMP1/xMx44dkZKSUmP5yJEjMXLkSKtljz/+uOX306dPx/Tp02t87uqxWQCs5sNyJAYrchqSJEGvU0OvU+MG/5rfTrmSscKM80W1fxuy+tfzl5ebKgUKSkwoKDHh2HWmpNCqFDUG3gd46GqMA/Nz10DFKSmIiFocBiuiWmhUCoR4uSDE6/pTUhSUmP7+NmRR2TXnByssq0B5hRlnLpbizMX6T0lR6zgwd+spKYiIqHlgsCJqBEmS4O2mgbebBh3rOSVF7uWB+FdfAasOYeeLjDZPSeHvoYHSpMT/DPsRqNdZwpi//u8rY75uWigVvA1JRORIDFZETaShU1LUdiuy+lVUXoFSUyVO55cCkJB+MOea7SokwNfdevD9tW5FckoKouajtvFCZH/2Os4MVkTNjC1TUpQYqwbjn7tYjPVbdiCsQzdcKL56tvy/p6SoDmTX46FV1ZgZv7YpKbxc1FDwKhiRQyiVVf/BMRqNcHGpe1gCNV5JSdV8idd7LNz1MFgRtWCuGhXCfVUI0WuQ6ytwV9+2tf6lUFFpRn6xsc7B+DWmpCivwMnzjZ+SojqccUoKItuoVCq4uroiLy8ParX6urOpm81mGI1GlJWVyT7zuhwa2n8hBEpKSpCbmwsvLy9LoG0oBisiJ6BSKqpmpNfr6tzumlNSFJUjz/LeAVNSWK6EcUoKomqSJCE4OBjp6enIyMi47vZCCJSWlsLFxcUpf4Ya238vL68GPQT6agxWRGTRkCkpavsGpGWZoWowvrHS3OgpKf5+SDenpCDnodFo0KFDBxiNxutuazKZsGXLFgwYMKDRt7Naosb0X61WN/pKVTUGKyJqEFumpLhUarpiJnzrKSmqb0vmFZbD0MApKaqDVmGOArkpGQj2crV6YLe7ln/VUculUCig09V9tRmoGpNVUVEBnU7nlMGqufSff9sQkUNJkgQvVw28XBs+JUX1bcnrT0mhwMZzR2q066pR1rjiVdtVMR83DaekIKJGYbAiomajvlNSVJoFLpYYrcaBZReUYNefR+HmG4wLxSbLlbBiYyVKjJXIuFCCjAsldbarVEjwddP8Pe6r+hak1eSsVd+Q1Kk5GJ+IamKwIqIWR3nFlBTVTCYT2hb/hbvu6mF1G6C4vKLGbcerJ2bNKyzDheKqq2C5l5dfj4dOZfUcSKuHc18xJszLVe2UA4mJnBWDFRG1am5aFdy0KkT4udW5XUWlGReKjX8Pwq/l0UTV34osrzCjsKwChWUVOJl3/Skpqh9BdOV0FFbfjvSo+r1GxcH4RC0dgxUREaqmpAjU6xCo1wHwvOZ21VNS/B28yqxmw78yiFVPSXHuUhnOXSoDcKnOGrxc1X/P/XWtKSncddC7cEoKouaKwYqIyAZXTknRPqBhU1L8fWuy9ikpjubUPSWFRqWoMQmrr6sa2TkSdEfyEOzligAPHXzdNVBzSgqiJsVgRUTkIA2ZksIqfBlqn5LCeM0pKZRYfnKf5d3VU1IE1DIzPqekILIv/iQREcms4VNSVA28zy0sR86lUvx5IhMKVy/kFZXXMSVF7a6ckuLKEMYpKYhsw2BFRNSCXGtKCpPJhLVrM3DXXX2hVqthNgvklxivmAvs78H31o8oKmvQlBQ1vgHJKSmIAMgcrLZs2YLZs2djz549yMrKwo8//oiRI0da1gshMG3aNHzxxRcoKCjAbbfdhs8++wwdOnQAACQnJ2PQoEG1tr1z507ccsstKCsrwxNPPIE9e/bg8OHDuPvuu7Fq1aoa2ycnJ2PSpEk4ePAgwsLCMGXKFDzyyCNW28ydOxezZ89GdnY2evTogY8//hi33nqrvQ4HEZHdKK6YkqJLcN3bXjklRa3jwK4xJcXB69TgoVXBX1/7lBRXXhXz5pQU1IrIGqyKi4vRo0cP/POf/8R9991XY31iYiI++ugjLF68GJGRkXj99dcRFxeHQ4cOQafTISYmBllZWVafef3117Fp0yb07t0bAFBZWQkXFxc899xzWLlyZa11pKen4x//+AeeeOIJLF26FJs2bcL48eMRHByMuLg4AMDy5csxadIkzJs3D3369MEHH3yAuLg4HDlyBAEBAXY+MkRETaehU1JcOQ6s1ikpyitQmNfwKSmuHgfGKSmoJZA1WA0bNgzDhg2rdZ0QAh988AGmTJmCe+65BwCwZMkSBAYGYtWqVRg9ejQ0Go3Vk6hNJhN++uknPPvss5b//bi5ueGzzz4DAGzbtg0FBQU19jVv3jxERkbivffeAwB06dIFv//+O95//31LsJozZw4ef/xxPProo5bPrFmzBgsWLMCrr75qnwNCRNSM2WNKiquvijVmSoqrB+P7uqpQYbZvn4ls1WzHWKWnpyM7OxuxsbGWZZ6enujTpw9SUlIwevToGp9ZvXo1Lly4YAk/9ZWSkmK1HwCIi4vDxIkTAQBGoxF79uzBa6+9ZlmvUCgQGxuLlJSUa7ZbXl6O8vK/Z3A2GAwAqgKgyWSyqcbrqW7P3u22FOy/c/cf4DFobv13UQLh3lqEe2sB6K+5XXmFGReKypFXZMT5wnLkXp6eIq/IePnXcuQVGnG+qBymSnHdKSlu9VdgWDM5Bk2tuZ0DTc2R/belzWYbrLKzswEAgYGBVssDAwMt6642f/58xMXFoU2bNjbvq7b9GAwGlJaW4uLFi6isrKx1m7/++uua7c6aNQszZsyosXzDhg1wda37WWgNlZSU5JB2Wwr237n7D/AYtPT+e15+tVcB8L78AiAEUFIBGEyAwShd/hUwmCTLr8cNEnafl/D9miR4auTrg9xa+jnQWI7of0lJ3V/quFKzDVa2OnPmDNavX48VK1bIXYrFa6+9hkmTJlneGwwGhIWFYejQodDrr/0/uIYwmUxISkrCkCFDrJ6T5izYf+fuP8Bj4Oz9B4BRn6dib+Yl5Lq3x5jYjnKX0+Sc/RxwZP+r7zjVR7MNVtVjp3JychAc/PdXWnJyctCzZ88a2y9cuBC+vr4YMWJEg/aVk5NjtSwnJwd6vR4uLi5QKpVQKpW1bnPlGK+rabVaaLXaGsvVarXDTnpHtt0SsP/O3X+Ax8CZ+z+ub1vszTyAFXuz8PzQLk4767wznwOAY/pvS3vN9qyLjIxEUFAQNm3aZFlmMBiQmpqK6Ohoq22FEFi4cCHi4+MbdDCjo6Ot9gNUXUqs3o9Go0FUVJTVNmazGZs2bapRCxERyWNo10Do1VXTQaw/WPuQESJHk/WKVVFREY4fP255n56ejrS0NPj4+KBt27aYOHEi3nzzTXTo0MEy3UJISIjVXFcAsHnzZqSnp2P8+PG17ufQoUMwGo3Iz89HYWEh0tLSAMBy5euJJ57AJ598gsmTJ+Of//wnNm/ejBUrVmDNmjWWNiZNmoSEhAT07t0bt956Kz744AMUFxfbPFCeiIgcQ6NSICZQYN0ZCUu2Z+Du7iFyl0ROSNZgtXv3bqsJPqvHIyUkJGDRokWYPHkyiouLMWHCBBQUFKBfv35Yt24ddDqdVTvz589HTEwMOnfuXOt+7rrrLmRkZFje9+rVC0DVlS6g6urYmjVr8MILL+DDDz9EmzZt8OWXX1qmWgCAUaNGIS8vD1OnTkV2djZ69uyJdevW1RjQTkRE8okJNGPjOSV2nsrHoXMGdA2x73hWouuRNVgNHDjQEm5qI0kSZs6ciZkzZ9bZzrJly+pcf+rUqXrVsm/fvjq3eeaZZ/DMM89cty0iIpKHpwaI6xqINX9mY0nKKfzn/u5yl0ROptmOsSIiImqIh/uGAQBWpZ1FQYlR5mrI2TBYERFRqxLV1gtdgvUoM5nx3e4zcpdDTobBioiIWhVJkpAQHQ4A+GpHBirN1x5yQmRvDFZERNTq3NMzFJ4uapzOL8FvR3PlLoecCIMVERG1Oi4aJUbdUjXWavH2jOtsTWQ/DFZERNQqPdwnHJIE/HY0D+nni+Uuh5wEgxUREbVKbX1dcUenAADAVym8akVNg8GKiIharfiYCADAd3syUVxeIW8x5BQYrIiIqNXq394PkX5uKCyrwI/7zspdDjkBBisiImq1FAoJ4/pWTb2wJOVUnU/7ILIHBisiImrV7o9qA1eNEkdzirDjZL7c5VArx2BFREStmqeLGvf2CgVQddWKyJEYrIiIqNWLj44AAGw4lINzBaXyFkOtGoMVERG1ep2CPNC3nQ8qzQLLUk/LXQ61YgxWRETkFBIuX7X6ZudplFdUylsMtVoMVkRE5BSGdA1EsKcOF4qNWHsgS+5yqJVisCIiIqegUirw8OWpF/j8QHIUBisiInIao24Jg0apQFpmAfZnFshdDrVCDFZEROQ0/Ny1uLt7MABgMadeIAdgsCIiIqdS/fzAX/Zn4UJRubzFUKvDYEVERE6lZ5gXerTxhLHSjG93ZcpdDrUyDFZEROR0qicMXbojAxWVZnmLoVaFwYqIiJzOP7oHw8dNg3OXyrDxcK7c5VArwmBFREROR6dWYvQtYQD4/ECyLwYrIiJySg/3DYdCArafuIBjOYVyl0OtBIMVERE5pRAvFwztGgQAWJLCCUPJPhisiIjIacXHVM3EvnLvGRjKTDJXQ60BgxURETmt6Ha+6BDgjhJjJVbuOSN3OdQKMFgREZHTkiTJMmHoVykZMJuFvAVRi8dgRURETu2+XqHw0Kpw8nwxfj9+Xu5yqIVjsCIiIqfmplXh/qg2ADj1AjUegxURETm9cdFVg9g3/ZWLzPwSmauhlozBioiInN4N/u7o38EPQgBf7+DUC9RwDFZEREQAHrk8iP3bXZkoNVbKWwy1WAxWREREAAZ2CkCYjwsulZrw8/5zcpdDLRSDFREREQClQsK4vlVjrRZtPwUhOPUC2Y7BioiI6LIHe4dBq1LgUJYBe09flLscaoEYrIiIiC7zctVgZM9QAMCi7RzETrZjsCIiIrpC9dQL/zuQhVxDmczVUEvDYEVERHSFG0M90TvcGxVmgWU7T8tdDrUwDFZERERXqX5+4NLU0zBWmOUthloUBisiIqKr3NktCP4eWuQVlmP9wWy5y6EWhMGKiIjoKhqVAmNvbQuAzw8k2zBYERER1WJsn7ZQKSTsOnURB89dkrscaiEYrIiIiGoRqNfhzhuDAABfpXDqBaofBisiIqJrqH5+4Kq0sygoMcpbDLUIDFZERETXEBXuja7BepSZzPhu9xm5y6EWgMGKiIjoGiRJQkJM1YShX+3IQKWZzw+kujFYERER1WFEj1B4uqhxOr8EyUdy5S6HmjkGKyIiojq4aJQYdUsYAGAxB7HTdTBYERERXcfDfcIhScCWo3k4mVckdznUjDFYERERXUdbX1fc0SkAQNVYK6JrYbAiIiKqh+rnB36/+wyKyyvkLYaaLQYrIiKieujf3g+Rfm4oLK/Aj/vOyl0ONVMMVkRERPWgUEgY17dq6oUlKacgBKdeoJoYrIiIiOrp/qg2cNUocTSnCDtO5stdDjVDDFZERET15Omixr29QgFUXbUiupqswWrLli0YPnw4QkJCIEkSVq1aZbVeCIGpU6ciODgYLi4uiI2NxbFjxyzrk5OTIUlSra9du3ZZtvvjjz/Qv39/6HQ6hIWFITExsUYtH3zwATp16gQXFxeEhYXhhRdeQFlZmdU2c+fORUREBHQ6Hfr06YOdO3fa94AQEVGzFx8dAQDYcCgH5wpK5S2Gmh1Zg1VxcTF69OiBuXPn1ro+MTERH330EebNm4fU1FS4ubkhLi7OEnhiYmKQlZVl9Ro/fjwiIyPRu3dvAIDBYMDQoUMRHh6OPXv2YPbs2Zg+fTo+//xzy36WLVuGV199FdOmTcPhw4cxf/58LF++HP/+978t2yxfvhyTJk3CtGnTsHfvXvTo0QNxcXHIzeUsvEREzqRTkAf6tvNBpVlgWeppucuhZkYl586HDRuGYcOG1bpOCIEPPvgAU6ZMwT333AMAWLJkCQIDA7Fq1SqMHj0aGo0GQUFBls+YTCb89NNPePbZZyFJEgBg6dKlMBqNWLBgATQaDbp164a0tDTMmTMHEyZMAABs374dt912G8aOHQsAiIiIwJgxY5Cammppe86cOXj88cfx6KOPAgDmzZuHNWvWYMGCBXj11Vdr7UN5eTnKy8st7w0Gg6VOk8nUoGN2LdXt2bvdloL9d+7+AzwGzt5/oGmPwUO3hmHHyXws25mBJwZEQKuSf2SNs58Djuy/LW3KGqzqkp6ejuzsbMTGxlqWeXp6ok+fPkhJScHo0aNrfGb16tW4cOGCJfwAQEpKCgYMGACNRmNZFhcXh3feeQcXL16Et7c3YmJi8PXXX2Pnzp249dZbcfLkSaxduxbjxo0DABiNRuzZswevvfaapQ2FQoHY2FikpKRcsw+zZs3CjBkzaizfsGEDXF1dbTsg9ZSUlOSQdlsK9t+5+w/wGDh7/4GmOQaVAvDSKJFfbMJ/lq7HLf7N5xuCzn4OOKL/JSUl9d622Qar7OxsAEBgYKDV8sDAQMu6q82fPx9xcXFo06aNVTuRkZE12qhe5+3tjbFjx+L8+fPo168fhBCoqKjAE088YbkVeP78eVRWVtZay19//XXNPrz22muYNGmS5b3BYEBYWBiGDh0KvV5/vUNgE5PJhKSkJAwZMgRqtdqubbcE7L9z9x/gMXD2/gNNfwzOuJ/EnI3H8UeZN6bd1dfh+7seZz8HHNn/6jtO9dFsg5Wtzpw5g/Xr12PFihU2fzY5ORlvv/02Pv30U/Tp0wfHjx/H888/jzfeeAOvv/56g2vSarXQarU1lqvVaoed9I5suyVg/527/wCPgbP3H2i6YzC2bwQ++fUk/jhjwMHsYvQM83L4PuvD2c8BR/Tflvbkvyl8DdVjp3JycqyW5+TkWI2rqrZw4UL4+vpixIgRNdqprY0r9/H6669j3LhxGD9+PG666Sbce++9ePvttzFr1iyYzWb4+flBqVTWuxYiImr9/Ny1uLt7MABOvUB/a7bBKjIyEkFBQdi0aZNlmcFgQGpqKqKjo622FUJg4cKFiI+Pr5Eqo6OjsWXLFquBZ0lJSejUqRO8vb0BVN07VSisD4VSqbS0rdFoEBUVZVWL2WzGpk2batRCRETOo/r5gb/sz8KFovK6NyanIGuwKioqQlpaGtLS0gBUDVhPS0vD6dOnIUkSJk6ciDfffBOrV6/GgQMHEB8fj5CQEIwcOdKqnc2bNyM9PR3jx4+vsY+xY8dCo9Hgsccew8GDB7F8+XJ8+OGHVmOfhg8fjs8++wzffvst0tPTkZSUhNdffx3Dhw+3BKxJkybhiy++wOLFi3H48GE8+eSTKC4uthooT0REzqVnmBd6tPGEsdKMb3dlyl0ONQOyjrHavXs3Bg0aZHlfHXYSEhKwaNEiTJ48GcXFxZgwYQIKCgrQr18/rFu3Djqdzqqd+fPnIyYmBp07d66xD09PT2zYsAFPP/00oqKi4Ofnh6lTp1qmWgCAKVOmQJIkTJkyBWfPnoW/vz+GDx+Ot956y7LNqFGjkJeXh6lTpyI7Oxs9e/bEunXragxoJyIi5xIfHYEXv9uPpTsy8K8B7aBSNtubQdQEZA1WAwcOrPMhlpIkYebMmZg5c2ad7SxbtqzO9d27d8fWrVuvuV6lUmHatGmYNm1ane0888wzeOaZZ+rchoiInMs/ugfjrbWHce5SGTYezsWdN3LsrTNjrCYiImoEnVqJ0beEAeAgdmKwIiIiarSH+oZDIQHbT1zAsZxCucshGTFYERERNVKolwuGdK0ac7skJUPmakhODFZERER2kBAdAQBYufcMDGXO+bw+YrAiIiKyi+gbfNEhwB0lxkqs3HNG7nJIJgxWREREdiBJkmXC0K9SMmA2N58HM1PTYbAiIiKyk/t6hcJDq8LJ88X4/fh5ucshGTBYERER2YmbVoX7o9oA4NQLzorBioiIyI7GRYcDADb9lYvM/BKZq6GmxmBFRERkRzf4u6N/Bz8IAXy9g1MvOBsGKyIiIjurnnrh212ZKDVWylsMNSkGKyIiIjsb1DkAbbxdcKnUhNX7z8pdDjUhBisiIiI7UyokjOtbNdZq8fYMCMGpF5wFgxUREZEDPNg7DFqVAoeyDNh7+qLc5VATYbAiIiJyAG83De7pGQKg6qoVOQcGKyIiIgeJvzyIfe2BLOQayuQthpoEgxUREZGD3Bjqid7h3qgwCyzbeVrucqgJMFgRERE5UPXzA5emnoaxwixvMeRwDFZEREQOdGe3IPh7aJFXWI71B7PlLoccjMGKiIjIgTQqBcbe2hYAnx/oDBisiIiIHGxsn7ZQKSTsOnURB89dkrsccqAGByuj0YgzZ87g9OnTVi8iIiKyFqjX4c4bgwAAX6Vw6oXWzOZgdezYMfTv3x8uLi4IDw9HZGQkIiMjERERgcjISEfUSERE1OIlXB7EvirtLApKjPIWQw6jsvUDjzzyCFQqFX755RcEBwdDkiRH1EVERNSq9A73RtdgPQ5lGfDd7jN4fEA7uUsiB7A5WKWlpWHPnj3o3LmzI+ohIiJqlSRJQkJMOF5ZeQBf7cjAP/tFQqngxYnWxuZbgV27dsX58+cdUQsREVGrNqJHKDxd1DidX4LkI7lyl0MOYHOweueddzB58mQkJyfjwoULMBgMVi8iIiKqnYtGiVG3hAEAFnMQe6tk863A2NhYAMDgwYOtlgshIEkSKisr7VMZERFRK/Rwn3B8sfUkthzNw8m8IrTzd5e7JLIjm4PVr7/+6og6iIiInEJbX1fc0SkAm/7KxVc7MjBteDe5SyI7silYmUwmzJw5E/PmzUOHDh0cVRMREVGrFh8TgU1/5eL73Wfw0tBOcNPafJ2Dmimbxlip1Wr88ccfjqqFiIjIKfRv74dIPzcUllfgx31n5S6H7MjmwesPP/ww5s+f74haiIiInIJCIWFc33AAVc8PFELIXBHZi83XHisqKrBgwQJs3LgRUVFRcHNzs1o/Z84cuxVHRETUWj3Quw3e3XAER3OKsONkPqJv8JW7JLIDm4PVn3/+iZtvvhkAcPToUat1nIWdiIiofvQ6Ne67ORRf7ziNxdtPMVi1EvxWIBERkUzioyPw9Y7TSDqcg3MFpQjxcpG7JGokm8dYERERkX10DPRAdDtfVJoFlqZywtDWwOYrVoMGDarzlt/mzZsbVRAREZEzSYgJR8rJC/hmZyaevaMDdGql3CVRI9gcrHr27Gn13mQyIS0tDX/++ScSEhLsVRcREZFTiO0SiGBPHbIulWHtgSzcd3MbuUuiRrA5WL3//vu1Lp8+fTqKiooaXRAREZEzUSkVeLhvOGavP4LFKRkMVi2c3cZYPfzww1iwYIG9miMiInIao24Jg0apwP7MAqRlFshdDjWC3YJVSkoKdDqdvZojIiJyGn7uWtzdPRhA1YSh1HLZfCvwvvvus3ovhEBWVhZ2796N119/3W6FEREROZOEmAj8sO8sftmfhf+7qwt83bVyl0QNYPMVK71eD09PT8vLx8cHAwcOxNq1azFt2jRH1EhERNTq9QjzQo8wLxgrzfh2V6bc5VAD2XzFatGiRQ4og4iIiBKiwzEpswBLd2TgXwPaQaXkdJMtjc1/Yu3atcOFCxdqLC8oKEC7du3sUhQREZEzuuumYPi6aXDuUhk2Hs6VuxxqAJuD1alTp1BZWVljeXl5Oc6ePWuXooiIiJyRTq3E6FvDAHAQe0tV71uBq1evtvx+/fr18PT0tLyvrKzEpk2bEBERYdfiiIiInM1DfcLxWfIJbD9xAUdzCtEx0EPuksgG9Q5WI0eOBABIklRjhnW1Wo2IiAi89957di2OiIjI2YR4uWBo1yCsO5iNJSmn8ObIm+QuiWxQ71uBZrMZZrMZbdu2RW5uruW92WxGeXk5jhw5grvvvtuRtRIRETmF+JhwAMAPe8/CUGaSuRqyhc1jrNLT0+Hn5wcAKCsrs3tBREREzi66nS86BLijxFiJlXvOyF0O2cDmYGU2m/HGG28gNDQU7u7uOHnyJADg9ddfx/z58+1eIBERkbORJAnxMREAgK9SMmA2C3kLonqzOVi9+eabWLRoERITE6HRaCzLb7zxRnz55Zd2LY6IiMhZ3dcrFB5aFU6eL8bvx8/LXQ7Vk83BasmSJfj888/x0EMPQalUWpb36NEDf/31l12LIyIiclZuWhUe6N0GAKdeaElsDlZnz55F+/btayw3m80wmTjAjoiIyF7G9a0axL7pr1xk5pfIXA3Vh83BqmvXrti6dWuN5d9//z169epll6KIiIgIaOfvjgEd/SEE8PWODLnLoXqw+VmBU6dORUJCAs6ePQuz2YwffvgBR44cwZIlS/DLL784okYiIiKnlRAdji1H8/DtrkxMjO0IF43y+h8i2dh8xeqee+7Bzz//jI0bN8LNzQ1Tp07F4cOH8fPPP2PIkCGOqJGIiMhpDewUgDAfF1wqNWH1fj46rrmzKVhVVFRg5syZiIyMRFJSEnJzc1FSUoLff/8dQ4cOtXnnW7ZswfDhwxESEgJJkrBq1Sqr9UIITJ06FcHBwXBxcUFsbCyOHTtmWZ+cnAxJkmp97dq1y7LdH3/8gf79+0On0yEsLAyJiYk1aikoKMDTTz+N4OBgaLVadOzYEWvXrrXaZu7cuYiIiIBOp0OfPn2wc+dOm/tMRERkC6VCsoy1Wrw9A0Jw6oXmzKZgpVKpkJiYiIqKCrvsvLi4GD169MDcuXNrXZ+YmIiPPvoI8+bNQ2pqKtzc3BAXF2eZmDQmJgZZWVlWr/HjxyMyMhK9e/cGABgMBgwdOhTh4eHYs2cPZs+ejenTp+Pzzz+37MdoNGLIkCE4deoUvv/+exw5cgRffPEFQkNDLdssX74ckyZNwrRp07B371706NEDcXFxyM3l08eJiMixHuwdBq1KgUNZBuzJuCh3OVQHm8dYDR48GL/99ptdHrg8bNgwDBs2rNZ1Qgh88MEHmDJlCu655x4AVVM9BAYGYtWqVRg9ejQ0Gg2CgoIsnzGZTPjpp5/w7LPPQpIkAMDSpUthNBqxYMECaDQadOvWDWlpaZgzZw4mTJgAAFiwYAHy8/Oxfft2qNVqAKjRvzlz5uDxxx/Ho48+CgCYN28e1qxZgwULFuDVV1+ttQ/l5eUoLy+3vDcYDJY67f0Nyur2nPWbmey/c/cf4DFw9v4DrfsYuKkljOgRjO/2nMXCbenoEVrzwcytuf/14cj+29KmJGy8pjhv3jzMmDEDDz30EKKiouDm5ma1fsSIEbY093chkoQff/zR8rDnkydP4oYbbsC+ffvQs2dPy3a33347evbsiQ8//LBGGytXrsSDDz6IjIwMtGlTNfdHfHw8DAaD1W3GX3/9FXfccQfy8/Ph7e2Nu+66Cz4+PnB1dcVPP/0Ef39/jB07Fq+88gqUSiWMRiNcXV3x/fffW+oDgISEBBQUFOCnn36qtU/Tp0/HjBkzaixftmwZXF1dbT9IRETktM4UA7P/UEEhCUy/uRKemut/huyjpKQEY8eOxaVLl6DX6+vc1uYrVk899RSAqis4V5MkCZWVlbY2Wavs7GwAQGBgoNXywMBAy7qrzZ8/H3FxcZZQVd1OZGRkjTaq13l7e+PkyZPYvHkzHnroIaxduxbHjx/HU089BZPJhGnTpuH8+fOorKystZa6JkV97bXXMGnSJMt7g8GAsLAwDB069Lp/MLYymUxISkrCkCFDLFfdnAn779z9B3gMnL3/gHMcg80FO7HndAHO6zthzB03WK1zhv7XxZH9r77jVB82Byuz2WzrR5rEmTNnsH79eqxYscLmz5rNZgQEBODzzz+HUqlEVFQUzp49i9mzZ2PatGkNrkmr1UKr1dZYrlarHXbSO7LtloD9d+7+AzwGzt5/oHUfg4TbIrHn9D58s/sMnhncERpVzaHSrbn/9eGI/tvSns3TLTSV6rFTOTk5VstzcnKsxlVVW7hwIXx9fWvcigwKCqq1jSv3ERwcjI4dO1o9oqdLly7Izs6G0WiEn58flEplvWshIiJyhDu7BcHfQ4u8wnKsP1j73RuSV7MNVpGRkQgKCsKmTZssywwGA1JTUxEdHW21rRACCxcuRHx8fI1UGR0djS1btlgNPEtKSkKnTp3g7e0NALjttttw/Phxq6txR48eRXBwMDQaDTQaDaKioqxqMZvN2LRpU41aiIiIHEWjUmDsrW0B8PmBzZWswaqoqAhpaWlIS0sDAKSnpyMtLQ2nT5+GJEmYOHEi3nzzTaxevRoHDhxAfHw8QkJCrAaQA8DmzZuRnp6O8ePH19jH2LFjodFo8Nhjj+HgwYNYvnw5PvzwQ6uxT08++STy8/Px/PPP4+jRo1izZg3efvttPP3005ZtJk2ahC+++AKLFy/G4cOH8eSTT6K4uNjyLUEiIqKm8FCftlApJOw6dREHz12Suxy6is1jrOxp9+7dGDRokOV9ddhJSEjAokWLMHnyZBQXF2PChAkoKChAv379sG7dOuh0Oqt25s+fj5iYGHTu3LnGPjw9PbFhwwY8/fTTiIqKgp+fH6ZOnWqZagEAwsLCsH79erzwwgvo3r07QkND8fzzz+OVV16xbDNq1Cjk5eVh6tSpyM7ORs+ePbFu3boaA9qJiIgcKUCvw7CbgvHz/nNYsj0D7zzQXe6S6AqyBquBAwfWOYOsJEmYOXMmZs6cWWc7y5Ytq3N99+7da31w9JWio6OxY8eOOrd55pln8Mwzz9S5DRERkaMlRIfj5/3nsCrtLF67qzO8XDn3QnPRoGBlNptx/Phx5Obm1viW4IABA+xSGBEREdUuKtwbXYP1OJRlwIrdmZgw4Ibrf4iahM3BaseOHRg7diwyMmo+r8ie81gRERFR7SRJQkJMOF5ZeQBf7cjAY/3ayV0SXWbz4PUnnngCvXv3xp9//on8/HxcvHjR8srPz3dEjURERHSVET1C4emiRmZ+KZKP8Lm1zYXNV6yOHTuG77//Hu3bt3dEPURERFQPLholRt0Shs+3nMTilAwMaO8jd0mEBlyx6tOnD44fP+6IWoiIiMgGD/cJhyQBW47mIf18sdzlEBpwxerZZ5/Fiy++iOzsbNx00001JuTs3p1f+yQiImoKbX1dcUenAGz6Kxdfp2YiSpK7IrI5WN1///0AgH/+85+WZZIkQQjBwetERERNLD4mApv+ysUP+86hWw+5qyGbg1V6eroj6iAiIqIG6N/eD5F+bkg/X4zdeRLuk7sgJ2dzsAoPD3dEHURERNQACoWEcX3DMfOXQ9iarahz4m1yvAbPvH7o0CGcPn0aRqPRavmIESMaXRQRERHV3wO92+DdDUeQXVqJ1PSL6N+Jj1uTi83B6uTJk7j33ntx4MABy9gqoGqcFQCOsSIiImpiep0a9/QIxje7zuCr1NMMVjKyebqF559/HpGRkcjNzYWrqysOHjyILVu2oHfv3khOTnZAiURERHQ9D/cJAwBsPJyLswWlMlfjvGwOVikpKZg5cyb8/PygUCigUCjQr18/zJo1C88995wjaiQiIqLr6BjogfZ6M8wCWJaaIXc5TsvmYFVZWQkPDw8AgJ+fH86dOwegalD7kSNH7FsdERER1duAoKrhOd/szESZiUNz5GBzsLrxxhuxf/9+AFWzsCcmJmLbtm2YOXMm2rXjQyCJiIjkcqOPQJBei/xiI9YeyJK7HKdkc7CaMmUKzGYzAGDmzJlIT09H//79sXbtWnz00Ud2L5CIiIjqRykBY2+tGmu1OIW3A+Vg87cC4+LiLL9v3749/vrrL+Tn58Pb29vyzUAiIiKSx4NRofjk15PYn1mAtMwC9Azzkrskp2LzFatqx48fx/r161FaWgofHz5Rm4iIqDnwddfi7u7BAIAlKafkLcYJ2RysLly4gMGDB6Njx4646667kJVVdQ/3sccew4svvmj3AomIiMg28TERAIBf9mfhfFG5vMU4GZuD1QsvvAC1Wo3Tp0/D1dXVsnzUqFFYt26dXYsjIiIi2/UM80KPNp4wVpqxfFem3OU4FZuD1YYNG/DOO++gTZs2Vss7dOiAjAwOlCMiImoO4qMjAABf78hARaVZ3mKciM3Bqri42OpKVbX8/HxotVq7FEVERESN84/uwfBx0yDrUhk2Hs6RuxynYXOw6t+/P5YsWWJ5L0kSzGYzEhMTMWjQILsWR0RERA2jUysx+pbLUy9s5x2lpmLzdAuJiYkYPHgwdu/eDaPRiMmTJ+PgwYPIz8/Htm3bHFEjERERNcBDfcMx77cTSDl5AUdzCtEx0EPuklq9Bs28fvToUfTr1w/33HMPiouLcd9992Hfvn244YYbHFEjERERNUColwuGdA0EwKkXmorNV6wAwNPTE//3f/9n71qIiIjIzhJiIrD+YA5+2HsWk+/sDL1OLXdJrVqDglVZWRn++OMP5ObmWh5vU23EiBF2KYyIiIgaL7qdLzoEuONYbhFW7jmDR2+LlLukVs3mYLVu3TrEx8fj/PnzNdZJkoTKSj5Nm4iIqLmQJAnxMRF4fdWf+ColAwnREVAo+Ag6R7F5jNWzzz6L//f//h+ysrJgNputXgxVREREzc99vULhoVXh5Pli/H685oURsh+bg1VOTg4mTZqEwMBAR9RDREREduamVeH+qKqJvTmI3bFsDlYPPPAAkpOTHVAKEREROcq46HAAwKa/cnH6QonM1bReNo+x+uSTT/D//t//w9atW3HTTTdBrbb+dsFzzz1nt+KIiIjIPm7wd0f/Dn7Yeuw8lqZm4LW7ushdUqtkc7D65ptvsGHDBuh0OiQnJ0OS/h4AJ0kSgxUREVEzlRAdga3HzuPbXZmYGNsRLhql3CW1OjbfCvy///s/zJgxA5cuXcKpU6eQnp5ueZ08edIRNRIREZEdDOocgDbeLrhUasLq/WflLqdVsjlYGY1GjBo1CgqFzR8lIiIiGSkVEsb1rRprtXh7BoQQMlfU+ticjhISErB8+XJH1EJEREQONuqWMGhVChzKMmBPxkW5y2l1bB5jVVlZicTERKxfvx7du3evMXh9zpw5diuOiIiI7MvLVYORPUOxfHcmFqdkoHeEj9wltSo2B6sDBw6gV69eAIA///zTat2VA9mJiIioeRoXHY7luzPxvwNZyP1HFwTodXKX1GrYHKx+/fVXR9RBRERETeTGUE/0DvfG7oyLWLbzNCbGdpS7pFaDI9CJiIicUHxMBABgaeppGCvM8hbTijBYEREROaE7uwXB30OLvMJyrDuYLXc5rQaDFRERkRPSqBQYe2tbAMCS7afkLaYVYbAiIiJyUg/1aQuVQsLujIv48+wluctpFRisiIiInFSAXodhNwUDAL5KyZC5mtaBwYqIiMiJJURXzcS+Ku0sCkqMMlfT8jFYERERObGocG90DdajvMKMFbsz5S6nxWOwIiIicmKSJCEhpuqq1Vc7MlBp5vMDG4PBioiIyMmN6BEKTxc1MvNLkXwkV+5yWjQGKyIiIifnolFi1C1hAIBFnHqhURisiIiICOP6hkOSgK3HzuNEXpHc5bRYDFZERESEMB9XDO4cAIBTLzQGgxUREREBAOKjIwAAK/ecQVF5hbzFtFAMVkRERAQA6NfeD+383FBYXoEf956Ru5wWicGKiIiIAAAKhYRxlycMXZKSASE49YKtGKyIiIjI4v6oNnDVKHEstwgpJy/IXU6Lw2BFREREFnqdGvfdHAoAWLKdg9htJWuw2rJlC4YPH46QkBBIkoRVq1ZZrRdCYOrUqQgODoaLiwtiY2Nx7Ngxy/rk5GRIklTra9euXZbt/vjjD/Tv3x86nQ5hYWFITEy8Zk3ffvstJEnCyJEjbaqFiIiotagexL7hUDbOFpTKW0wLI2uwKi4uRo8ePTB37txa1ycmJuKjjz7CvHnzkJqaCjc3N8TFxaGsrAwAEBMTg6ysLKvX+PHjERkZid69ewMADAYDhg4divDwcOzZswezZ8/G9OnT8fnnn9fY36lTp/DSSy+hf//+NtdCRETUWnQM9EB0O1+YBbAslVetbCFrsBo2bBjefPNN3HvvvTXWCSHwwQcfYMqUKbjnnnvQvXt3LFmyBOfOnbNc2dJoNAgKCrK8fH198dNPP+HRRx+FJEkAgKVLl8JoNGLBggXo1q0bRo8ejeeeew5z5syx2l9lZSUeeughzJgxA+3atbO5FiIiotYkISYCAPDNzkyUmSrlLaYFUcldwLWkp6cjOzsbsbGxlmWenp7o06cPUlJSMHr06BqfWb16NS5cuIBHH33UsiwlJQUDBgyARqOxLIuLi8M777yDixcvwtvbGwAwc+ZMBAQE4LHHHsPWrVsbXQsAlJeXo7y83PLeYDAAAEwmE0wmky2H47qq27N3uy0F++/c/Qd4DJy9/wCPgb37f3t7bwR76pB1qQyr953Bvb1C7NKuozjyz9+WNpttsMrOzgYABAYGWi0PDAy0rLva/PnzERcXhzZt2li1ExkZWaON6nXe3t74/fffMX/+fKSlpdmtFgCYNWsWZsyYUWP5hg0b4Orqes3PNUZSUpJD2m0p2H/n7j/AY+Ds/Qd4DOzZ/yhPCb9cUuKT9QegzUqzW7uO5Ig//5KSknpv22yDla3OnDmD9evXY8WKFTZ9rrCwEOPGjcMXX3wBPz8/u9b02muvYdKkSZb3BoMBYWFhGDp0KPR6vV33ZTKZkJSUhCFDhkCtVtu17ZaA/Xfu/gM8Bs7ef4DHwBH971NsxPrZv+F0MRDa/Tb0aONpl3YdwZF//tV3nOqj2QaroKAgAEBOTg6Cg4Mty3NyctCzZ88a2y9cuBC+vr4YMWJEjXZycnKsllW/DwoKwokTJ3Dq1CkMHz7cst5sNgMAVCoVjhw5YnMt1bRaLbRabY3larXaYT/0jmy7JWD/nbv/AI+Bs/cf4DGwZ/+DvNQY3iMEP+w9i2U7z6B3pH0vQDiCI/78bWmv2c5jFRkZiaCgIGzatMmyzGAwIDU1FdHR0VbbCiGwcOFCxMfH1+h8dHQ0tmzZYnV/NCkpCZ06dYK3tzc6d+6MAwcOIC0tzfIaMWIEBg0ahLS0NISFhdlUCxERUWuScHnqhV/+yML5ovK6NyZ5g1VRUZElzABVg8TT0tJw+vRpSJKEiRMn4s0338Tq1atx4MABxMfHIyQkpMYcU5s3b0Z6ejrGjx9fYx9jx46FRqPBY489hoMHD2L58uX48MMPLbfodDodbrzxRquXl5cXPDw8cOONN0Kj0dhUCxERUWvSI8wLPcK8YKw0Y/muTLnLafZkvRW4e/duDBo0yPK+OuwkJCRg0aJFmDx5MoqLizFhwgQUFBSgX79+WLduHXQ6nVU78+fPR0xMDDp37lxjH56entiwYQOefvppREVFwc/PD1OnTsWECRNsqrW+tRAREbU2CdHhmJRZgK93ZOBfA9pBpWy2N7xkJ2uwGjhwYJ0PeJQkCTNnzsTMmTPrbGfZsmV1ru/evXuNKRTqsmjRogbXQkRE1NrcdVMw3lpzGFmXyrDxcA7uvDH4+h9yUoycREREVCedWonRt4YBABbz+YF1YrAiIiKi63qoTzgUEpBy8gKO5hTKXU6zxWBFRERE1xXi5YKhXaumH1qSckreYpoxBisiIiKql+rnB/6w9ywMZc756KDrYbAiIiKieunbzgcdA91RYqzE97vPyF1Os8RgRURERPUiSRLiL08Y+tWODJjN1/5mv7NisCIiIqJ6u7dXKDy0KqSfL8bW4+flLqfZYbAiIiKienPTqvBA7zYAgMXbT8lbTDPEYEVEREQ2Gdc3HADw65FcnL5QInM1zQuDFREREdmknb87BnT0hxDA16mcMPRKDFZERERks4ToqqtWy3dlotRYKXM1zQeDFREREdlsYKcAhPm44FKpCav3n5W7nGaDwYqIiIhsplRIlrFWi7dnQAhOvQAwWBEREVEDPdg7DFqVAoeyDNiTcVHucpoFBisiIiJqEC9XDUb2DAUALOLUCwAYrIiIiKgR4mOqbgeu+zMbOYYymauRH4MVERERNVi3EE/cEuGNCrPAstTTcpcjOwYrIiIiapTq5wcu23kaxgqzvMXIjMGKiIiIGiWuWxACPLTIKyzHuoPZcpcjKwYrIiIiahSNSoGxfdoCAJY4+SB2BisiIiJqtLG3toVKIWF3xkX8efaS3OXIhsGKiIiIGi1Ar8Owm4IBAF+lOO/zAxmsiIiIyC6qnx+4Ku0sCkqMMlcjDwYrIiIisouocG90DdajvMKMFbsz5S5HFgxWREREZBeSJCHh8oShX+3IQKXZ+Z4fyGBFREREdjOiRyg8XdTIzC9F8pFcuctpcgxWREREZDcuGiVG3RIGwDmfH8hgRURERHb1cJ9wSBKw9dh5nMgrkrucJsVgRURERHbV1tcVd3QKAOB8Uy8wWBEREZHdxcdEAAC+33MGReUV8hbThBisiIiIyO76t/dDpJ8bisor8OO+s3KX02QYrIiIiMjuFAoJ4/pWTb2wZPspCOEcUy8wWBEREZFDPNC7DVw1ShzLLULKyQtyl9MkGKyIiIjIIfQ6Ne67ORQAsGS7cwxiZ7AiIiIih4mPjgAAbDiUjbMFpfIW0wQYrIiIiMhhOgZ6ILqdL8wCWJba+q9aMVgRERGRQ1U/P/CbnZkoM1XKXI1jMVgRERGRQ8V2CUSwpw75xUasPZAldzkOxWBFREREDqVSKvDw5akXFrfy5wcyWBEREZHDjbolDBqlAvvPXEJaZoHc5TgMgxURERE5nJ+7Fnd3DwZQNWFoa8VgRURERE2i+vmBv/yRhfNF5fIW4yAMVkRERNQkeoZ5oUcbTxgrzVi+K1PuchyCwYqIiIiaTMLlq1Zf78hARaVZ3mIcgMGKiIiImsxdNwXD102DrEtl2Hg4R+5y7I7BioiIiJqMTq3E6FvDAACLW+HzAxmsiIiIqEk91CccCglIOXkBR3MK5S7HrhisiIiIqEmFeLlgaNcgAMCSlFPyFmNnDFZERETU5OIvPz/wh71nYSgzyVyN/TBYERERUZOLbueLDgHuKDFWYuWeM3KXYzcMVkRERNTkJEmyTBi6JCUDZrOQtyA7YbAiIiIiWdzXKxQeWhXSzxdj6/HzcpdjFwxWREREJAs3rQr3R7UB0HqeH8hgRURERLKJj64axL75SC5OXyiRuZrGY7AiIiIi2bTzd8eAjv4QAvg6teVPGMpgRURERLJKuHzVavmuTJQaK2WupnEYrIiIiEhWAzsFIMzHBZdKTVi9/6zc5TSKrMFqy5YtGD58OEJCQiBJElatWmW1XgiBqVOnIjg4GC4uLoiNjcWxY8cs65OTkyFJUq2vXbt2Wbb7448/0L9/f+h0OoSFhSExMdFqP1988QX69+8Pb29veHt7IzY2Fjt37rSpFiIiImoYpULCuL5VV60Wb8+AEC136gVZg1VxcTF69OiBuXPn1ro+MTERH330EebNm4fU1FS4ubkhLi4OZWVlAICYmBhkZWVZvcaPH4/IyEj07t0bAGAwGDB06FCEh4djz549mD17NqZPn47PP//csp/k5GSMGTMGv/76K1JSUhAWFoahQ4fi7Nmz9a6FiIiIGu7B3mHQqhQ4lGXAnoyLcpfTYCo5dz5s2DAMGzas1nVCCHzwwQeYMmUK7rnnHgDAkiVLEBgYiFWrVmH06NHQaDQICgqyfMZkMuGnn37Cs88+C0mSAABLly6F0WjEggULoNFo0K1bN6SlpWHOnDmYMGGCZZsrffnll1i5ciU2bdqE+Pj4etVCREREDeflqsHInqFYvjsTi1My0DvCR+6SGkTWYFWX9PR0ZGdnIzY21rLM09MTffr0QUpKSq1hZvXq1bhw4QIeffRRy7KUlBQMGDAAGo3GsiwuLg7vvPMOLl68CG9v7xrtlJSUwGQywcfHp8G1AEB5eTnKy8st7w0GA4CqAGgy2fe5SNXt2bvdloL9d+7+AzwGzt5/gMegNfR/7K1Vwep/B7JwZmh7BOp19f6sI/tvS5vNNlhlZ2cDAAIDA62WBwYGWtZdbf78+YiLi0ObNm2s2omMjKzRRvW62oLVK6+8gpCQEEuQakgtADBr1izMmDGjxvINGzbA1dX1mp9rjKSkJIe021Kw/87df4DHwNn7D/AYtPT+t/NQ4mQh8OY3yRgWZrb5847of0lJ/efXarbBylZnzpzB+vXrsWLFika185///AfffvstkpOTodPVPynX5rXXXsOkSZMs7w0Gg2X8ll6vb1TbVzOZTEhKSsKQIUOgVqvt2nZLwP47d/8BHgNn7z/AY9Ba+i/CsjFxxR/YU6DDe48NgEZVv+Hgjux/9R2n+mi2wap67FROTg6Cg4Mty3NyctCzZ88a2y9cuBC+vr4YMWJEjXZycnKsllW/v3J8FgC8++67+M9//oONGzeie/fuDa6lmlarhVarrbFcrVY77KR3ZNstAfvv3P0HeAycvf8Aj0FL7/9d3UPx9v+OILewHJuOXsCIHiE2fd4R/belvWY7j1VkZCSCgoKwadMmyzKDwYDU1FRER0dbbSuEwMKFCxEfH1+j89HR0diyZYvV/dGkpCR06tTJ6jZgYmIi3njjDaxbt87yjcKG1EJEREQNp1EpMLZPWwAt8/mBsgaroqIipKWlIS0tDUDVIPG0tDScPn0akiRh4sSJePPNN7F69WocOHAA8fHxCAkJwciRI63a2bx5M9LT0zF+/Pga+xg7diw0Gg0ee+wxHDx4EMuXL8eHH35odYvunXfeweuvv44FCxYgIiIC2dnZyM7ORlFREQDYVAsRERE1zthb20KlkLA74yL+PHtJ7nJsIuutwN27d2PQoEGW99VhJyEhAYsWLcLkyZNRXFyMCRMmoKCgAP369cO6detqjH2aP38+YmJi0Llz5xr78PT0xIYNG/D0008jKioKfn5+mDp1qmWqBQD47LPPYDQa8cADD1h9dtq0aZg+fToA1LsWIiIiapwAvQ7DbgrGz/vP4auUDLzzQPfrf6iZkDVYDRw4sM7ZVSVJwsyZMzFz5sw621m2bFmd67t3746tW7dec/2pU6fq/LwttRAREVHjJUSH4+f957Aq7Sxeu6szvFw11/9QM9Bsx1gRERGR84oK90a3ED3KK8xYsTtT7nLqjcGKiIiImh1JkpAQHQEAWJKSgUpzy3h+IIMVERERNUsjeobAy1WNMxdL8etfuXKXUy8MVkRERNQs6dRKjOodBgBYnHJK3mLqicGKiIiImq2H+4ZDkoCtx87jRF6R3OVcF4MVERERNVthPq4Y3DkAAPBVSobM1VwfgxURERE1a/GXB7Gv3HMGReUV8hZzHQxWRERE1Kz1a++Hdn5uKCyvwI/7zspdTp0YrIiIiKhZUygkxEeHA6h6fmBdk4vLjcGKiIiImr37o9rATaPEsdwipJy8IHc518RgRURERM2eh06N+25uAwBYsr35DmJnsCIiIqIWofp24IZD2ThbUCpzNbVjsCIiIqIWoUOgB2Ju8IVZAEt3NM+rVgxWRERE1GJUT73w7a5MlJkq5S2mFgxWRERE1GLEdglAiKcO+cVGrPkjS+5yamCwIiIiohZDpVTgob6Xp15ohs8PZLAiIiKiFmX0LWHQqBTYf+YS0jIL5C7HCoMVERERtSi+7loM7x4CoGrC0OaEwYqIiIhanISYqtuBv/yRhfNF5TJX8zcGKyIiImpxurfxQs8wLxgrzVi+K1PuciwYrIiIiKhFqr5q9fWODFRUmmWupgqDFREREbVId90UDF83DbIulWHjX3lylwOAwYqIiIhaKK1KiTG3tgUAfJ16WuZqqjBYERERUYs1tk9bKBUSUtMv4lyJ3NUwWBEREVELFuLlgqFdAwEAv2fLH2vkr4CIiIioEaqfH7grT4Kh1CRrLQxWRERE1KL1beeDjgHuMJol/JB2TtZaGKyIiIioRZMkCQ/3DYNKEigokfeKlUrWvRMRERHZwcgeIVCeO4AHB7eXtQ5esSIiIqIWz0WjhLta7ioYrIiIiIjshsGKiIiIyE4YrIiIiIjshMGKiIiIyE4YrIiIiIjshMGKiIiIyE4YrIiIiIjshMGKiIiIyE4YrIiIiIjshMGKiIiIyE4YrIiIiIjshMGKiIiIyE4YrIiIiIjsRCV3Ac5ECAEAMBgMdm/bZDKhpKQEBoMBanUzeLx3E2P/nbv/AI+Bs/cf4DFg/x3X/+p/t6v/Ha8Lg1UTKiwsBACEhYXJXAkRERHZqrCwEJ6ennVuI4n6xC+yC7PZjHPnzsHDwwOSJNm1bYPBgLCwMGRmZkKv19u17ZaA/Xfu/gM8Bs7ef4DHgP13XP+FECgsLERISAgUirpHUfGKVRNSKBRo06aNQ/eh1+ud8geqGvvv3P0HeAycvf8AjwH775j+X+9KVTUOXiciIiKyEwYrIiIiIjthsGoltFotpk2bBq1WK3cpsmD/nbv/AI+Bs/cf4DFg/5tH/zl4nYiIiMhOeMWKiIiIyE4YrIiIiIjshMGKiIiIyE4YrIiIiIjshMFKJtOnT4ckSVavzp07W22TkpKCO+64A25ubtDr9RgwYABKS0sBAMnJyTU+X/3atWsXAODUqVO1rt+xY4fVfr777jt07twZOp0ON910E9auXdsq+g9UzZb77rvvomPHjtBqtQgNDcVbb71ltZ/k5GTcfPPN0Gq1aN++PRYtWuTw/gNNcwxq24ckSXBzc7PaT2s+B9avX4++ffvCw8MD/v7+uP/++3Hq1Cmr/chxDjRV/1esWIGePXvC1dUV4eHhmD17do1aWurPAAAcPXoU99xzD/z8/KDX69GvXz/8+uuvVm2cPn0a//jHP+Dq6oqAgAC8/PLLqKiosNqmJZ4DQP36/9xzzyEqKgparRY9e/astZY//vgD/fv3h06nQ1hYGBITE+3e39o0xTHYv38/xowZg7CwMLi4uKBLly748MMPa9Rit3NAkCymTZsmunXrJrKysiyvvLw8y/rt27cLvV4vZs2aJf7880/x119/ieXLl4uysjIhhBDl5eVWn83KyhLjx48XkZGRwmw2CyGESE9PFwDExo0brbYzGo2W/Wzbtk0olUqRmJgoDh06JKZMmSLUarU4cOBAi++/EEI8++yzolOnTuKnn34SJ0+eFLt37xYbNmywrD958qRwdXUVkyZNEocOHRIff/yxUCqVYt26dQ7tf1Mdg8LCwhrbdO3aVSQkJFj205rPgZMnTwqtVitee+01cfz4cbFnzx4xYMAA0atXL8t+5DoHmqL/a9euFSqVSnz22WfixIkT4pdffhHBwcHi448/lr3/9jgGQgjRoUMHcdddd4n9+/eLo0ePiqeeekq4urqKrKwsIYQQFRUV4sYbbxSxsbFi3759Yu3atcLPz0+89tprsh+Dpui/EFV/D37yySdi3LhxokePHjXquHTpkggMDBQPPfSQ+PPPP8U333wjXFxcxH//+1+H9l+IpjkG8+fPF88995xITk4WJ06cEF999ZVwcXFx2M8Bg5VMpk2bVusJXq1Pnz5iypQp9W7PaDQKf39/MXPmTMuy6mC1b9++a37uwQcfFP/4xz9q7Ptf//pXvffdEE3R/0OHDgmVSiX++uuva35u8uTJolu3blbLRo0aJeLi4uq974ZqimNwtbS0NAFAbNmyxbKsNZ8D3333nVCpVKKystKybPXq1UKSJMt/MOQ6B5qi/2PGjBEPPPCA1XYfffSRaNOmjSV8teSfgby8vBrns8FgEABEUlKSEKIqXCoUCpGdnW3Z5rPPPhN6vV6Ul5cLIVruOVCf/tdnf59++qnw9va2HA8hhHjllVdEp06d6teRRmjqY1DtqaeeEoMGDbK8t+c5wFuBMjp27BhCQkLQrl07PPTQQzh9+jQAIDc3F6mpqQgICEBMTAwCAwNx++234/fff79mW6tXr8aFCxfw6KOP1lg3YsQIBAQEoF+/fli9erXVupSUFMTGxloti4uLQ0pKih16WDdH9//nn39Gu3bt8MsvvyAyMhIREREYP3488vPzLdvI2X+g6c6Bal9++SU6duyI/v37W5a15nMgKioKCoUCCxcuRGVlJS5duoSvvvoKsbGxUKvVAFp3/8vLy6HT6ay2c3FxwZkzZ5CRkQGgZf8M+Pr6olOnTliyZAmKi4tRUVGB//73vwgICEBUVJSlfzfddBMCAwOt+mcwGHDw4EHLNi3xHKhP/+sjJSUFAwYMgEajsSyLi4vDkSNHcPHiRft19hrkOAaXLl2Cj4+P5b1dzwGboxjZxdq1a8WKFSvE/v37xbp160R0dLRo27atMBgMIiUlRQAQPj4+YsGCBWLv3r1i4sSJQqPRiKNHj9ba3rBhw8SwYcOsluXl5Yn33ntP7NixQ+zcuVO88sorQpIk8dNPP1m2UavVYtmyZVafmzt3rggICLB/p6/QFP3/17/+JbRarejTp4/YsmWL+PXXX0XPnj2t/pfSoUMH8fbbb1t9bs2aNQKAKCkpsX/Hr9AUx+BKpaWlwtvbW7zzzjtWy1vzOSCEEMnJySIgIEAolUoBQERHR4uLFy9a1st1DjRF///73/8KV1dXsXHjRlFZWSmOHDkiOnfuLACI7du3CyFa/s9AZmamiIqKEpIkCaVSKYKDg8XevXst6x9//HExdOhQq/0WFxcLAGLt2rVCiJZ9Dlyv/1e61tWhIUOGiAkTJlgtO3jwoAAgDh06ZNc+X62pj4EQVcMfVCqVWL9+vWWZPc8BBqtm4uLFi0Kv14svv/xSbNu2TQCwGgMghBA33XSTePXVV2t8NjMzUygUCvH9999fdz/jxo0T/fr1s7yX6x/Vqzmi/48//rgAII4cOWJZtmfPHgHAcntQzn9Uruboc2DZsmVCpVJZ3RIRonWfA1lZWaJDhw7i5ZdfFnv37hW//fabuP3228XgwYMtt8KayzngiP6bzWYxefJkodPphFKpFN7e3mL69OkCgNixY4cQovn0Xwjbj4HZbBYjRowQw4YNE7///rvYs2ePePLJJ0VoaKg4d+6cEKJ5B6urOaL/V2qOwepqjj4GBw4cEH5+fuKNN96wWm7Pc4C3ApsJLy8vdOzYEcePH0dwcDAAoGvXrlbbdOnSxXKJ9EoLFy6Er68vRowYcd399OnTB8ePH7e8DwoKQk5OjtU2OTk5CAoKakg3GswR/Q8ODoZKpULHjh2t2gBgaeda/dfr9XBxcWl8x2zg6HPgyy+/xN133211SwRo3efA3Llz4enpicTERPTq1QsDBgzA119/jU2bNiE1NRVA8zkHHNF/SZLwzjvvoKioCBkZGcjOzsatt94KAGjXrh2A5tN/wPZjsHnzZvzyyy/49ttvcdttt+Hmm2/Gp59+ChcXFyxevBjAtftXva6ubZr7OVCf/tdHfY5RU3HkMTh06BAGDx6MCRMmYMqUKVbr7HkOMFg1E0VFRThx4gSCg4MRERGBkJAQHDlyxGqbo0ePIjw83GqZEAILFy5EfHy8ZcxIXdLS0iwnKwBER0dj06ZNVtskJSUhOjq6Eb2xnSP6f9ttt6GiogInTpywagOApZ3m0n/AsedAeno6fv31Vzz22GM11jWXY+CI/peUlEChsP5rTqlUAgDMZjOA1t3/akqlEqGhodBoNPjmm28QHR0Nf39/AM2n/4Dtx6CkpAQAavwZKxQKqz/fAwcOIDc317I+KSkJer3e8g92czkGjuh/fURHR2PLli0wmUyWZUlJSejUqRO8vb0b2p0GcdQxOHjwIAYNGoSEhIQaU+4Adj4HbLq+RXbz4osviuTkZJGeni62bdsmYmNjhZ+fn8jNzRVCCPH+++8LvV4vvvvuO3Hs2DExZcoUodPpxPHjx63a2bhxowAgDh8+XGMfixYtEsuWLROHDx8Whw8fFm+99ZZQKBRiwYIFlm2q7zW/++674vDhw2LatGlN8lX7puh/ZWWluPnmm8WAAQPE3r17xe7du0WfPn3EkCFDLNtUf8X25ZdfFocPHxZz585tsq+aN8UxqDZlyhQREhIiKioqaqxrzefApk2bhCRJYsaMGeLo0aNiz549Ii4uToSHh1su78t1DjRF//Py8sRnn30mDh8+LPbt2yeee+45odPpRGpqqmWblvwzkJeXJ3x9fcV9990n0tLSxJEjR8RLL70k1Gq1SEtLE0L8Pd3C0KFDRVpamli3bp3w9/evdbqFlnYO1Kf/Qghx7NgxsW/fPvGvf/1LdOzYUezbt0/s27fP8i3AgoICERgYKMaNGyf+/PNP8e233wpXV9cmmW6hKY7BgQMHhL+/v3j44YetpnWo3ocQ9j0HGKxkMmrUKBEcHCw0Go0IDQ0Vo0aNqvEX5qxZs0SbNm2Eq6uriI6OFlu3bq3RzpgxY0RMTEyt+1i0aJHo0qWLcHV1FXq9Xtx6663iu+++q7HdihUrRMeOHYVGoxHdunUTa9assU8n69AU/RdCiLNnz4r77rtPuLu7i8DAQPHII4+ICxcuWG1TPahdo9GIdu3aiYULF9qlj9fTVMegsrJStGnTRvz73/++5jat+Rz45ptvRK9evYSbm5vw9/cXI0aMqBFC5DgHmqL/eXl5om/fvsLNzU24urqKwYMHW8ZWXakl/wzs2rVLDB06VPj4+AgPDw/Rt29fy9ipaqdOnRLDhg0TLi4uws/PT7z44ovCZDJZbdNSz4H69P/2228XAGq80tPTLdvs379f9OvXT2i1WhEaGir+85//OKzfV2qKYzBt2rRa+x8eHm7Vjr3OAUkIIWy/zkVEREREV+MYKyIiIiI7YbAiIiIishMGKyIiIiI7YbAiIiIishMGKyIiIiI7YbAiIiIishMGKyIiIiI7YbAiIiIishMGKyIiIiI7YbAiImoESZKwatUqucsgomaCwYqI6BqMRqPcJRBRC8NgRUQt0i+//AIvLy9UVlYCANLS0iBJEl599VXLNuPHj8fDDz9seb9y5Up069YNWq0WEREReO+996zajIiIwBtvvIH4+Hjo9XpMmDABRqMRzzzzDIKDg6HT6RAeHo5Zs2ZZtgeAe++9F5IkWd7X5syZMxgzZgx8fHzg5uaG3r17IzU1FQBw4sQJ3HPPPQgMDIS7uztuueUWbNy40erzn376KTp06ACdTofAwEA88MADlnVmsxmzZs1CZGQkXFxc0KNHD3z//fe2H1QiajSV3AUQETVE//79UVhYiH379qF379747bff4Ofnh+TkZMs2v/32G1555RUAwJ49e/Dggw9i+vTpGDVqFLZv346nnnoKvr6+eOSRRyyfeffddzF16lRMmzYNAPDRRx9h9erVWLFiBdq2bYvMzExkZmYCAHbt2oWAgAAsXLgQd955J5RKZa21FhUV4fbbb0doaChWr16NoKAg7N27F2az2bL+rrvuwltvvQWtVoslS5Zg+PDhOHLkCNq2bYvdu3fjueeew1dffYWYmBjk5+dj69atlvZnzZqFr7/+GvPmzUOHDh2wZcsWPPzww/D398ftt99uz8NORNcjiIhaqJtvvlnMnj1bCCHEyJEjxVtvvSU0Go0oLCwUZ86cEQDE0aNHhRBCjB07VgwZMsTq8y+//LLo2rWr5X14eLgYOXKk1TbPPvusuOOOO4TZbK61BgDixx9/rLPO//73v8LDw0NcuHCh3n3r1q2b+Pjjj4UQQqxcuVLo9XphMBhqbFdWViZcXV3F9u3brZY/9thjYsyYMfXeHxHZB28FElGLdfvttyM5ORlCCGzduhX33XcfunTpgt9//x2//fYbQkJC0KFDBwDA4cOHcdttt1l9/rbbbsOxY8cstxMBoHfv3lbbPPLII0hLS0OnTp3w3HPPYcOGDTbXmZaWhl69esHHx6fW9UVFRXjppZfQpUsXeHl5wd3dHYcPH8bp06cBAEOGDEF4eDjatWuHcePGYenSpSgpKQEAHD9+HCUlJRgyZAjc3d0tryVLluDEiRM210pEjcNbgUTUYg0cOBALFizA/v37oVar0blzZwwcOBDJycm4ePFig26Dubm5Wb2/+eabkZ6ejv/973/YuHEjHnzwQcTGxto0hsnFxaXO9S+99BKSkpLw7rvvon379nBxccEDDzxgGTzv4eGBvXv3Ijk5GRs2bMDUqVMxffp07Nq1C0VFRQCANWvWIDQ01KpdrVZb7xqJyD4YrIioxaoeZ/X+++9bQtTAgQPxn//8BxcvXsSLL75o2bZLly7Ytm2b1ee3bduGjh07XnNsVDW9Xo9Ro0Zh1KhReOCBB3DnnXciPz8fPj4+UKvVVle8atO9e3d8+eWXls9cbdu2bXjkkUdw7733Aqi6gnXq1CmrbVQqFWJjYxEbG4tp06bBy8sLmzdvxpAhQ6DVanH69GmOpyJqBhisiKjF8vb2Rvfu3bF06VJ88sknAIABAwbgwQcfhMlksgoaL774Im655Ra88cYbGDVqFFJSUvDJJ5/g008/rXMfc+bMQXBwMHr16gWFQoHvvvsOQUFB8PLyAlD1zcBNmzbhtttug1arhbe3d402xowZg7fffhsjR47ErFmzEBwcjH379iEkJATR0dHo0KEDfvjhBwwfPhySJOH111+3DGwHqr4BefLkSQwYMADe3t5Yu3YtzGYzOnXqBA8PD7z00kt44YUXYDab0a9fP1y6dAnbtm2DXq9HQkKCHY40EdWb3IO8iIga4/nnnxcAxOHDhy3LevToIYKCgmps+/3334uuXbsKtVot2rZtaxn4Xi08PFy8//77Vss+//xz0bNnT+Hm5ib0er0YPHiw2Lt3r2X96tWrRfv27YVKpRLh4eHXrPPUqVPi/vvvF3q9Xri6uorevXuL1NRUIYQQ6enpYtCgQcLFxUWEhYWJTz75RNx+++3i+eefF0IIsXXrVnH77bcLb29v4eLiIrp37y6WL19uadtsNosPPvhAdOrUSajVauHv7y/i4uLEb7/9Vt/DSER2IgkhhNzhjoiIiKg14LcCiYiIiOyEwYqIiIjIThisiIiIiOyEwYqIiIjIThisiIiIiOyEwYqIiIjIThisiIiIiOyEwYqIiIjIThisiIiIiOyEwYqIiIjIThisiIiIiOzk/wMx41s5iOlDxQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df = pd.DataFrame()\n", "\n", "for wc in np.linspace(56750, 56820, 50):\n", " m = farmer(crops, yields)\n", "\n", " m.eval(\"param wc;\")\n", " m.param[\"wc\"] = wc\n", " m.eval(\"s.t. min_profit{s in SCENARIOS}: scenario_profit[s] >= wc;\")\n", " # maximize mean profit\n", " m.eval(\n", " \"maximize objective: (sum{s in SCENARIOS} scenario_profit[s]) / card(SCENARIOS);\"\n", " )\n", " m.get_output(\"solve;\")\n", "\n", " df.loc[wc, \"return\"] = m.obj[\"objective\"].value()\n", " df.loc[wc, \"wc\"] = wc\n", "\n", "df.plot(x=\"wc\", y=\"return\", xlabel=\"worst case\", ylabel=\"mean return\", grid=True)" ] }, { "cell_type": "markdown", "id": "5e93ba71-3df5-472b-b9bf-a54c4be68237", "metadata": {}, "source": [ "What we see is that the worst case can be raised from 48,820 to 56,800 (a difference of 7,980 euro) at a cost of reducing the expected from from 108,390 to 107,100 (a reduction of 1,290 euro). This improvement in worst-case performance may be worth the reduction in mean profit." ] }, { "cell_type": "code", "execution_count": 14, "id": "54d3eb2d-d32b-4b27-bb95-26709f217d6a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Objective = 107100.00\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
profit
scenario
average117500.0
good147000.0
poor56800.0
\n", "
" ], "text/plain": [ " profit\n", "scenario \n", "average 117500.0\n", "good 147000.0\n", "poor 56800.0" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
acres
crop
beets300.0
corn100.0
wheat100.0
\n", "
" ], "text/plain": [ " acres\n", "crop \n", "beets 300.0\n", "corn 100.0\n", "wheat 100.0" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
growbuysellexcessexpenserevenueprofit
scenariocrop
averagebeets6000.006000.00.078000.0216000.0138000.0
corn300.0060.00.023000.09000.0-14000.0
wheat250.0050.00.015000.08500.0-6500.0
goodbeets7200.006000.01200.078000.0228000.0150000.0
corn360.00120.00.023000.018000.0-5000.0
wheat300.00100.00.015000.017000.02000.0
poorbeets4800.004800.00.078000.0172800.094800.0
corn240.000.00.023000.00.0-23000.0
wheat200.000.00.015000.00.0-15000.0
\n", "
" ], "text/plain": [ " grow buy sell excess expense revenue profit\n", "scenario crop \n", "average beets 6000.0 0 6000.0 0.0 78000.0 216000.0 138000.0\n", " corn 300.0 0 60.0 0.0 23000.0 9000.0 -14000.0\n", " wheat 250.0 0 50.0 0.0 15000.0 8500.0 -6500.0\n", "good beets 7200.0 0 6000.0 1200.0 78000.0 228000.0 150000.0\n", " corn 360.0 0 120.0 0.0 23000.0 18000.0 -5000.0\n", " wheat 300.0 0 100.0 0.0 15000.0 17000.0 2000.0\n", "poor beets 4800.0 0 4800.0 0.0 78000.0 172800.0 94800.0\n", " corn 240.0 0 0.0 0.0 23000.0 0.0 -23000.0\n", " wheat 200.0 0 0.0 0.0 15000.0 0.0 -15000.0" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "wc = 56800\n", "\n", "m = farmer(crops, yields)\n", "\n", "m.eval(\"param wc;\")\n", "m.param[\"wc\"] = wc\n", "m.eval(\"s.t. min_profit{s in SCENARIOS}: scenario_profit[s] >= wc;\")\n", "\n", "# maximize mean profit\n", "m.eval(\n", " \"maximize objective: (sum{s in SCENARIOS} scenario_profit[s]) / card(SCENARIOS);\"\n", ")\n", "m.get_output(\"solve;\")\n", "\n", "farm_report(m)" ] }, { "cell_type": "markdown", "id": "13e4cbb7-9248-424e-a825-4af1f4f89977", "metadata": {}, "source": [ "## Bibliographic Notes\n", "\n", "The Farmer's Problem is well known in the stochastic programming community. Examples of treatments include the following web resources:\n", "\n", "- http://user.engineering.uiowa.edu/~dbricker/stacks_pdf1/slpwr_farmer.pdf\n", "- https://www.math.uh.edu/~rohop/Spring_15/Chapter1.pdf" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.6" }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }