Farmer: Zero Value of the Stochastic Solution

import sys
from pathlib import Path
if "google.colab" in sys.modules:
    from urllib.request import urlretrieve
    base = "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/"
    for name in ["helper.py", "farmer.py"]:
        urlretrieve(base + name, name)
    import helper
    helper.easy_install()
else:
    sys.path.insert(0, str(Path.cwd().parent))
    import helper
import farmer as f
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display
helper.set_plotting_style()

from farmer import *

# Exercise 8 data, p. 19. Purchase prices are 40% above selling prices.
EX8_PLANT_COST = {"WHEAT": 180.0, "CORN": 280.0, "BEETS": 310.0}
EX8_SELL_PRICE = {
    "WHEAT": 300.0,
    "CORN": 170.0,
    "BEETS_FAVORABLE": 41.0,
    "BEETS_UNFAVORABLE": 11.0,
}


def price_scenario(wheat_price, yield_multiplier=1.0):
    """One scenario of Exercise 8: a wheat selling price, optionally with a yield."""
    sell = dict(EX8_SELL_PRICE, WHEAT=wheat_price)
    return {
        "yield": {c: MEAN_YIELD[c] * yield_multiplier for c in CROPS},
        "sell": sell,
        "buy": {c: 1.4 * sell[c] for c in FEED_CROPS},
    }


# Wait-and-see: what the farmer would do knowing next year's wheat price.
wait_and_see = {}
for price in (300.0, 220.0):
    m = solve_farmer(build_farmer([price_scenario(price)], plant_cost=EX8_PLANT_COST))
    wait_and_see[price] = pyo.value(m.expected_profit)
    p = plan_of(m) + 0.0  # normalize any -0.0 the LP returns
    print(
        f"perfect forecast, wheat at ${price:.0f}/T:"
        f"  acres {p[0]:5.1f} / {p[1]:5.1f} / {p[2]:5.1f}"
        f"   profit ${wait_and_see[price]:,.0f}"
    )

# Here-and-now: one plan, both prices, recourse after the price is revealed.
price_sp = solve_farmer(
    build_farmer(
        [price_scenario(300.0), price_scenario(220.0)], plant_cost=EX8_PLANT_COST
    )
)
price_plan = plan_of(price_sp)
rp_price = pyo.value(price_sp.expected_profit)
ws_price = np.mean(list(wait_and_see.values()))
print(
    f"\nstochastic plan: acres {price_plan[0]:5.1f} / {price_plan[1]:5.1f}"
    f" / {price_plan[2]:5.1f}   E[profit] ${rp_price:,.0f}"
)
print(f"WS   = ${ws_price:,.0f}")
print(f"RP   = ${rp_price:,.0f}")
print(f"EVPI = ${ws_price - rp_price:,.0f}   <- Exercise 8's question")

# EEV: plan at the mean price, then live with both.
mean_price_model = solve_farmer(
    build_farmer([price_scenario(260.0)], plant_cost=EX8_PLANT_COST)
)
eev, _ = grade_plan(
    plan_of(mean_price_model),
    [price_scenario(300.0), price_scenario(220.0)],
    plant_cost=EX8_PLANT_COST,
)
print(f"EEV  = ${eev:,.0f}")
print(f"VSS  = ${rp_price - eev:,.0f}")

def rows_of(record):
    rows=[]
    for s in record["scenarios"]:
        rows.append({"Scenario":s["id"],"Probability":s["probability"],
            **{"Buy "+k:v for k,v in s["purchases_tons"].items()},
            **{"Sell "+k:v for k,v in s["sales_tons"].items()},"Profit ($)":s["profit"]})
    display(pd.Series(record["acreage"],name="Acres").to_frame().round(3))
    display(pd.DataFrame(rows).round(2))
prices=[price_scenario(300),price_scenario(220)]
for price in [300,220]:
    print("Perfect information: wheat selling price",price)
    rows_of(solution_record(solve_farmer(build_farmer([price_scenario(price)],plant_cost=EX8_PLANT_COST))))
print("Shared stochastic plan")
rows_of(solution_record(price_sp))
print("Expected-value plan evaluated")
rows_of(evaluate_solution(plan_of(mean_price_model),prices,plant_cost=EX8_PLANT_COST))
assert np.isclose(rp_price-eev,0)
assert ws_price > rp_price