import%20marimo%0A%0A__generated_with%20%3D%20%220.23.16%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Time-of-Use%20Pricing%20for%20Residential%20Demand%20Response%0A%0A%20%20%20%20A%20utility%20would%20like%20to%20reduce%20the%20daily%20demand%20peak%2C%20but%20customers%20retain%0A%20%20%20%20control%20of%20when%20they%20consume%20flexible%20energy.%20This%20creates%20a%20bilevel%0A%20%20%20%20problem%3A%20the%20**upper%20problem**%20chooses%20a%2024-hour%20price%20vector%20for%20the%0A%20%20%20%20utility%2C%20while%20the%20**lower%20problem**%20models%20how%20a%20representative%20customer%0A%20%20%20%20reschedules%20consumption%20while%20balancing%20cost%20against%20discomfort.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20from%20pathlib%20import%20Path%0A%0A%20%20%20%20import%20cvxpy%20as%20cp%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20from%20blvpy%20import%20BilevelProblem%2C%20LowerProblem%0A%0A%20%20%20%20plt.style.use(Path(__file__).resolve().parents%5B1%5D%20%2F%20%22_shared%22%20%2F%20%22zhlatex.mplstyle%22)%0A%20%20%20%20return%20BilevelProblem%2C%20LowerProblem%2C%20Path%2C%20cp%2C%20mo%2C%20np%2C%20plt%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Customer%20response%20and%20utility%20objective%0A%0A%20%20%20%20Let%20%24p_t%24%20be%20the%20price%2C%20%24%5Cell_t%24%20the%20customer's%20load%2C%20and%0A%20%20%20%20%24%5Cbar%5Cell_t%24%20the%20preferred%20load%20in%20hour%20%24t%24.%20For%20a%20published%20price%20vector%2C%0A%20%20%20%20the%20customer%20solves%0A%0A%20%20%20%20%5C%5B%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%5Cmathop%7B%5Cmathrm%7Bminimize%7D%7D_%7B%5Cell%7D%0A%20%20%20%20%20%20%26%20p%5ET%5Cell%2B%5Cgamma%5ClVert%5Cell-%5Cbar%5Cell%5CrVert_2%5E2%20%5C%5C%0A%20%20%20%20%5Cmathop%7B%5Cmathrm%7Bsubject%5C%20to%7D%7D%0A%20%20%20%20%20%20%26%20%5Cmathbf%7B1%7D%5ET%5Cell%3DE%2C%5C%5C%0A%20%20%20%20%20%20%26%20%5Cell_%7B%5Cmin%7D%5Cleq%5Cell_t%5Cleq%5Cell_%7B%5Cmax%7D.%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%5C%5D%0A%0A%20%20%20%20The%20energy%20equality%20models%20shifting%20rather%20than%20curtailment.%20The%20quadratic%0A%20%20%20%20term%20represents%20inconvenience%20from%20departing%20from%20the%20preferred%20schedule%0A%20%20%20%20and%20makes%20the%20lower%20response%20unique.%0A%0A%20%20%20%20The%20upper%20problem%20introduces%20a%20peak%20epigraph%20%24z%24%20and%20anticipates%20the%20lower%0A%20%20%20%20response%3A%0A%0A%20%20%20%20%5C%5B%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%5Cmathop%7B%5Cmathrm%7Bminimize%7D%7D_%7Bp%2Cz%2C%5Cell%7D%0A%20%20%20%20%20%20%26%20z%2B%5Ceta%5ClVert%20p%5CrVert_2%5E2%2B%5Ckappa%5ClVert%20Dp%5CrVert_2%5E2%5C%5C%0A%20%20%20%20%5Cmathop%7B%5Cmathrm%7Bsubject%5C%20to%7D%7D%0A%20%20%20%20%20%20%26%20%5Cell_t%5Cleq%20z%2C%5Cquad%200%5Cleq%20z%5Cleq%204%2C%5C%5C%0A%20%20%20%20%20%20%26%200%5Cleq%20p_t%5Cleq%200.8%2C%5C%5C%0A%20%20%20%20%20%20%26%20%5Cell%5Ctext%7B%20solves%20the%20customer%20problem%20for%20%7Dp.%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%5C%5D%0A%0A%20%20%20%20The%20magnitude%20penalty%20discourages%20unnecessarily%20high%20prices%2C%20while%20the%0A%20%20%20%20difference%20operator%20%24D%24%20discourages%20abrupt%20hour-to-hour%20tariff%20changes.%0A%20%20%20%20The%20loose%20upper%20bound%20on%20%24z%24%20supplies%20a%20practical%20scale%20for%20this%20synthetic%0A%20%20%20%20example%3B%20the%20optimized%20peak%20lies%20well%20below%20it.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Synthetic%20daily%20profile%0A%0A%20%20%20%20The%20preferred%20schedule%20combines%20a%20modest%20morning%20rise%20with%20a%20larger%0A%20%20%20%20evening%20peak.%20All%2024%20hourly%20decisions%20remain%20in%20the%20model.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20hours%20%3D%20np.arange(24%2C%20dtype%3Dfloat)%0A%20%20%20%20preferred_load%20%3D%20(%0A%20%20%20%20%20%20%20%200.75%0A%20%20%20%20%20%20%20%20%2B%200.75%20*%20np.exp(-0.5%20*%20np.square((hours%20-%208.0)%20%2F%201.8))%0A%20%20%20%20%20%20%20%20%2B%201.65%20*%20np.exp(-0.5%20*%20np.square((hours%20-%2019.0)%20%2F%202.1))%0A%20%20%20%20)%0A%20%20%20%20daily_energy%20%3D%20float(np.sum(preferred_load))%0A%20%20%20%20minimum_hourly_load%20%3D%200.45%0A%20%20%20%20maximum_hourly_load%20%3D%203.2%0A%20%20%20%20discomfort_weight%20%3D%200.45%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20daily_energy%2C%0A%20%20%20%20%20%20%20%20discomfort_weight%2C%0A%20%20%20%20%20%20%20%20hours%2C%0A%20%20%20%20%20%20%20%20maximum_hourly_load%2C%0A%20%20%20%20%20%20%20%20minimum_hourly_load%2C%0A%20%20%20%20%20%20%20%20preferred_load%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Specify%20and%20solve%20the%20bilevel%20model%0A%0A%20%20%20%20%60price%60%20is%20the%20utility's%20vector%20decision%20and%20is%20listed%20as%20a%20lower-level%0A%20%20%20%20parameter.%20The%20%60.value%60%20assignments%20provide%20a%20neutral%20solving%20initial%20point.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20BilevelProblem%2C%0A%20%20%20%20LowerProblem%2C%0A%20%20%20%20cp%2C%0A%20%20%20%20daily_energy%2C%0A%20%20%20%20discomfort_weight%2C%0A%20%20%20%20maximum_hourly_load%2C%0A%20%20%20%20minimum_hourly_load%2C%0A%20%20%20%20np%2C%0A%20%20%20%20preferred_load%2C%0A)%3A%0A%20%20%20%20price%20%3D%20cp.Variable(24%2C%20name%3D%22hourly_price%22)%0A%20%20%20%20peak_load%20%3D%20cp.Variable(name%3D%22peak_load%22)%0A%20%20%20%20load%20%3D%20cp.Variable(24%2C%20name%3D%22hourly_load%22)%0A%20%20%20%20price.value%20%3D%20np.full(24%2C%200.4)%0A%20%20%20%20peak_load.value%20%3D%202.0%0A%0A%20%20%20%20lower_problem%20%3D%20LowerProblem(%0A%20%20%20%20%20%20%20%20cp.Minimize(price%20%40%20load%20%2B%20discomfort_weight%20*%20cp.sum_squares(load%20-%20preferred_load))%2C%0A%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20cp.sum(load)%20%3D%3D%20daily_energy%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20load%20%3E%3D%20minimum_hourly_load%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20load%20%3C%3D%20maximum_hourly_load%2C%0A%20%20%20%20%20%20%20%20%5D%2C%0A%20%20%20%20%20%20%20%20parameters%3D%5Bprice%5D%2C%0A%20%20%20%20)%0A%20%20%20%20upper_objective%20%3D%20peak_load%20%2B%200.025%20*%20cp.sum_squares(price)%20%2B%200.12%20*%20cp.sum_squares(price%5B1%3A%5D%20-%20price%5B%3A-1%5D)%0A%20%20%20%20problem%20%3D%20BilevelProblem(%0A%20%20%20%20%20%20%20%20cp.Minimize(upper_objective)%2C%0A%20%20%20%20%20%20%20%20lower_problem%2C%0A%20%20%20%20%20%20%20%20upper_constraints%3D%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20price%20%3E%3D%200.0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20price%20%3C%3D%200.8%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20peak_load%20%3E%3D%200.0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20peak_load%20%3C%3D%204.0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20load%20%3C%3D%20peak_load%2C%0A%20%20%20%20%20%20%20%20%5D%2C%0A%20%20%20%20)%0A%20%20%20%20return%20load%2C%20peak_load%2C%20price%2C%20problem%0A%0A%0A%40app.cell%0Adef%20_(problem)%3A%0A%20%20%20%20epsilon_target%20%3D%201e-5%0A%20%20%20%20result%20%3D%20problem.solve(%0A%20%20%20%20%20%20%20%20epsilon_initial%3D1e-2%2C%0A%20%20%20%20%20%20%20%20epsilon_target%3Depsilon_target%2C%0A%20%20%20%20%20%20%20%20verbose%3DFalse%2C%0A%20%20%20%20)%0A%20%20%20%20diagnostics%20%3D%20problem.gap_diagnostics(result)%0A%20%20%20%20return%20diagnostics%2C%20result%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Flat-price%20reference%20response%0A%0A%20%20%20%20With%20a%20zero%2C%20flat%20price%20the%20energy%20payment%20does%20not%20encourage%20shifting.%0A%20%20%20%20We%20solve%20that%20fixed-tariff%20customer%20problem%20independently%20with%20CVXPY%20to%0A%20%20%20%20obtain%20a%20reference%20peak%20and%20utility%20objective.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20cp%2C%0A%20%20%20%20daily_energy%2C%0A%20%20%20%20discomfort_weight%2C%0A%20%20%20%20maximum_hourly_load%2C%0A%20%20%20%20minimum_hourly_load%2C%0A%20%20%20%20np%2C%0A%20%20%20%20preferred_load%2C%0A)%3A%0A%20%20%20%20baseline_load_variable%20%3D%20cp.Variable(24%2C%20name%3D%22baseline_hourly_load%22)%0A%20%20%20%20baseline_price%20%3D%20np.zeros(24)%0A%20%20%20%20baseline_problem%20%3D%20cp.Problem(%0A%20%20%20%20%20%20%20%20cp.Minimize(%0A%20%20%20%20%20%20%20%20%20%20%20%20baseline_price%20%40%20baseline_load_variable%0A%20%20%20%20%20%20%20%20%20%20%20%20%2B%20discomfort_weight%20*%20cp.sum_squares(baseline_load_variable%20-%20preferred_load)%0A%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20cp.sum(baseline_load_variable)%20%3D%3D%20daily_energy%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20baseline_load_variable%20%3E%3D%20minimum_hourly_load%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20baseline_load_variable%20%3C%3D%20maximum_hourly_load%2C%0A%20%20%20%20%20%20%20%20%5D%2C%0A%20%20%20%20)%0A%20%20%20%20baseline_problem.solve(solver%3Dcp.CLARABEL)%0A%20%20%20%20baseline_load%20%3D%20np.asarray(baseline_load_variable.value%2C%20dtype%3Dfloat)%0A%20%20%20%20baseline_peak%20%3D%20float(np.max(baseline_load))%0A%20%20%20%20return%20baseline_load%2C%20baseline_peak%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(%0A%20%20%20%20baseline_peak%2C%0A%20%20%20%20daily_energy%2C%0A%20%20%20%20diagnostics%2C%0A%20%20%20%20load%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20np%2C%0A%20%20%20%20peak_load%2C%0A%20%20%20%20price%2C%0A%20%20%20%20result%2C%0A)%3A%0A%20%20%20%20assert%20result.succeeded%2C%20result.message%0A%0A%20%20%20%20optimized_load%20%3D%20np.asarray(load.value%2C%20dtype%3Dfloat)%0A%20%20%20%20optimized_price%20%3D%20np.asarray(price.value%2C%20dtype%3Dfloat)%0A%20%20%20%20optimized_peak%20%3D%20float(np.max(optimized_load))%0A%20%20%20%20assert%20abs(float(np.sum(optimized_load))%20-%20daily_energy)%20%3C%3D%201e-5%2C%20(%0A%20%20%20%20%20%20%20%20%22The%20optimized%20load%20does%20not%20preserve%20total%20daily%20energy.%22%0A%20%20%20%20)%0A%20%20%20%20assert%20optimized_peak%20%3C%20baseline_peak%20-%200.2%2C%20%22The%20optimized%20tariff%20did%20not%20materially%20reduce%20peak%20demand.%22%0A%0A%20%20%20%20mo.md(f%22%22%22%0A%20%20%20%20%23%23%20Result%0A%0A%20%20%20%20%7C%20quantity%20%7C%20value%20%7C%0A%20%20%20%20%7C%20---%20%7C%20---%3A%20%7C%0A%20%20%20%20%7C%20status%20%7C%20%60%7Bresult.status%7D%60%20%7C%0A%20%20%20%20%7C%20flat-price%20peak%20load%20%7C%20%7Bbaseline_peak%3A.4f%7D%20%7C%0A%20%20%20%20%7C%20optimized%20peak%20load%20%7C%20%7Boptimized_peak%3A.4f%7D%20%7C%0A%20%20%20%20%7C%20peak%20epigraph%20%7C%20%7Bfloat(peak_load.value)%3A.4f%7D%20%7C%0A%20%20%20%20%7C%20utility%20objective%20%7C%20%7Bfloat(result.objective)%3A.6f%7D%20%7C%0A%20%20%20%20%7C%20final%20epsilon%20%7C%20%7Bresult.final_epsilon%3A.3e%7D%20%7C%0A%20%20%20%20%7C%20maximum%20lifted%20residual%20%7C%20%7Bresult.residuals.max_violation%3A.3e%7D%20%7C%0A%20%20%20%20%7C%20complementarity%20%7C%20%7Bresult.complementarity%3A.3e%7D%20%7C%0A%20%20%20%20%7C%20signed%20source%20gap%20%7C%20%7Bdiagnostics.source_gap%3A.3e%7D%20%7C%0A%0A%20%20%20%20The%20tariff%20raises%20prices%20around%20the%20evening%20peak%20and%20shifts%20energy%20into%0A%20%20%20%20lower-demand%20hours%20without%20changing%20total%20daily%20consumption.%20The%20optimized%0A%20%20%20%20peak%20is%20materially%20below%20the%20flat-price%20reference%2C%20even%20after%20penalizing%0A%20%20%20%20tariff%20magnitude%20and%20roughness.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%20optimized_load%2C%20optimized_price%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Load%20shifting%20and%20the%20price%20signal%0A%0A%20%20%20%20The%20upper%20panel%20compares%20the%20customer's%20preferred%20schedule%20with%20its%0A%20%20%20%20optimized%20response.%20The%20lower%20panel%20shows%20the%20time-of-use%20price%20that%0A%20%20%20%20induces%20the%20shift.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20Path%2C%0A%20%20%20%20baseline_load%2C%0A%20%20%20%20hours%2C%0A%20%20%20%20optimized_load%2C%0A%20%20%20%20optimized_price%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20preferred_load%2C%0A)%3A%0A%20%20%20%20figure_dir%20%3D%20Path(__file__).resolve().parent%20%2F%20%22figures%22%0A%20%20%20%20figure_dir.mkdir(parents%3DTrue%2C%20exist_ok%3DTrue)%0A%0A%20%20%20%20fig%2C%20axes%20%3D%20plt.subplots(%0A%20%20%20%20%20%20%20%202%2C%0A%20%20%20%20%20%20%20%201%2C%0A%20%20%20%20%20%20%20%20figsize%3D(8.0%2C%206.0)%2C%0A%20%20%20%20%20%20%20%20sharex%3DTrue%2C%0A%20%20%20%20%20%20%20%20gridspec_kw%3D%7B%22height_ratios%22%3A%20%5B2.0%2C%201.0%5D%7D%2C%0A%20%20%20%20)%0A%20%20%20%20axes%5B0%5D.plot(hours%2C%20preferred_load%2C%20linestyle%3D%22--%22%2C%20linewidth%3D2.0%2C%20label%3Dr%22%24%5Cbar%7B%5Cell%7D%24%22)%0A%20%20%20%20axes%5B0%5D.plot(hours%2C%20baseline_load%2C%20linestyle%3D%22%3A%22%2C%20linewidth%3D1.5%2C%20label%3D%22Flat-price%20response%22)%0A%20%20%20%20axes%5B0%5D.plot(hours%2C%20optimized_load%2C%20linewidth%3D2.2%2C%20label%3D%22Optimized-price%20response%22)%0A%20%20%20%20axes%5B0%5D.set(ylabel%3Dr%22%24%5Cell%24%22)%0A%20%20%20%20axes%5B0%5D.grid(alpha%3D0.25)%0A%20%20%20%20axes%5B0%5D.legend(frameon%3DFalse)%0A%0A%20%20%20%20axes%5B1%5D.bar(hours%2C%20optimized_price%2C%20width%3D0.85%2C%20color%3D%22%23e41a1c%22%2C%20alpha%3D0.85)%0A%20%20%20%20axes%5B1%5D.set(%0A%20%20%20%20%20%20%20%20xlabel%3D%22%24t%24%22%2C%0A%20%20%20%20%20%20%20%20ylabel%3D%22%24p%24%22%2C%0A%20%20%20%20%20%20%20%20xticks%3Drange(0%2C%2024%2C%203)%2C%0A%20%20%20%20)%0A%20%20%20%20axes%5B1%5D.grid(axis%3D%22y%22%2C%20alpha%3D0.25)%0A%0A%20%20%20%20fig.tight_layout()%0A%20%20%20%20figure_path%20%3D%20figure_dir%20%2F%20%22demand_response_pricing.pdf%22%0A%20%20%20%20fig.savefig(figure_path%2C%20bbox_inches%3D%22tight%22)%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A
7d6b68a6992ad84444bd8622891b273d