PROJECT 105 · PYTHON SOURCE

Electricity demand planning

Study-specific code. Shared modules, dependency versions, and reproduction instructions are included in all project files.

from pathlib import Path
import sys
sys.path.insert(0, str(Path(__file__).resolve().parents[2]))
"""Complete-hour energy accounting and calibrated next-day forecasts."""
from capstones.common import *
from capstones.fetch import SOURCES
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error
from threadpoolctl import threadpool_limits
import duckdb

SLUG='105-electricity-demand-planning'

def run():
    p=project_dir(SLUG);u,n=SOURCES[4];archive=download(u,n);provenance([(archive,u)],p)
    with zipfile.ZipFile(archive) as z:
        data=pd.read_csv(z.open('household_power_consumption.txt'),sep=';',na_values=['?'],low_memory=False)
    data['timestamp']=pd.to_datetime(data.Date+' '+data.Time,format='%d/%m/%Y %H:%M:%S')
    raw=len(data);missing=int(data.Global_active_power.isna().sum())
    con=duckdb.connect();con.register('minute_readings',data)
    sql="""-- Each minute-average kW contributes kW / 60 kWh.
-- Complete-hour targets require 60 present minute readings, not a partial sum.
SELECT date_trunc('hour',timestamp) AS timestamp,
count(*) AS calendar_minutes, count(Global_active_power) AS measured_minutes,
CASE WHEN count(Global_active_power)=60 THEN sum(Global_active_power)/60 END AS kwh,
sum(Global_active_power)/60 AS observed_kwh,
sum(Sub_metering_1)/1000 AS kitchen_kwh,
sum(Sub_metering_2)/1000 AS laundry_kwh,
sum(Sub_metering_3)/1000 AS water_heat_ac_kwh
FROM minute_readings GROUP BY 1 ORDER BY 1;
"""
    (p/'sql/hourly_mart.sql').write_text(sql,encoding='utf-8');hourly=con.execute(sql).df();con.close()
    frame=hourly.set_index('timestamp').reindex(pd.date_range(hourly.timestamp.min(),hourly.timestamp.max(),freq='h'))
    frame.index.name='timestamp'
    frame['hour']=frame.index.hour;frame['weekday']=frame.index.dayofweek;frame['month']=frame.index.month
    features=['hour','weekday','month']
    for lag in [24,48,168,336]:
        frame[f'lag_{lag}']=frame.kwh.shift(lag);features.append(f'lag_{lag}')
    frame['prior_week_same_hour']=pd.concat([frame.kwh.shift(24*day) for day in range(1,8)],axis=1).mean(axis=1)
    features.append('prior_week_same_hour')
    eligible=frame[(frame.index>='2007-01-01')&frame.kwh.notna()].copy()
    train=eligible[eligible.index<'2009-01-01'];val=eligible[(eligible.index>='2009-01-01')&(eligible.index<'2009-07-01')]
    calibration=eligible[(eligible.index>='2009-07-01')&(eligible.index<'2010-01-01')];test=eligible[eligible.index>='2010-01-01'].copy()
    median=float(train.kwh.median())
    def fallback(f): return f.lag_168.fillna(f.lag_24).fillna(f.prior_week_same_hour).fillna(median).to_numpy()
    baseline_candidates={'Previous week':lambda f:fallback(f),'Previous day':lambda f:f.lag_24.fillna(f.prior_week_same_hour).fillna(median).to_numpy(),'Prior seven days same hour':lambda f:f.prior_week_same_hour.fillna(median).to_numpy()}
    base_name=min(baseline_candidates,key=lambda name:mean_absolute_error(val.kwh,baseline_candidates[name](val)))
    configs=[];models=[]
    with threadpool_limits(limits=3):
        for leaves in [7,15]:
            model=HistGradientBoostingRegressor(loss='squared_error',max_leaf_nodes=leaves,max_iter=140,learning_rate=.06,l2_regularization=8,early_stopping=False,random_state=42)
            model.fit(train[features],train.kwh)
            configs.append({'leaves':leaves,'validation_mae_kwh':mean_absolute_error(val.kwh,np.maximum(0,model.predict(val[features])))})
            models.append(model)
        selected=models[int(np.argmin([c['validation_mae_kwh'] for c in configs]))]
        cal_prediction=np.maximum(0,selected.predict(calibration[features]));test['model']=np.maximum(0,selected.predict(test[features]))
    residual=np.abs(calibration.kwh.to_numpy()-cal_prediction)
    rank=min(len(residual)-1,int(np.ceil((len(residual)+1)*.9))-1);radius=float(np.sort(residual)[rank])
    test['baseline']=baseline_candidates[base_name](test);test['lower']=np.maximum(0,test.model-radius);test['upper']=test.model+radius
    test['covered']=((test.kwh>=test.lower)&(test.kwh<=test.upper)).astype(int)
    test['day']=test.index.strftime('%Y-%m-%d');test['period']=test.index.strftime('%Y-%m');test['absolute_error']=(test.model-test.kwh).abs();test['baseline_error']=(test.baseline-test.kwh).abs()
    results=pd.DataFrame([{'model':name,'test_hours':len(test),'mae_kwh':mean_absolute_error(test.kwh,pred),'rmse_kwh':np.sqrt(mean_squared_error(test.kwh,pred))} for name,pred in [(base_name,test.baseline),('Gradient boosting',test.model)]])
    save_csv(results,p/'outputs/results.csv');save_csv(hourly,p/'data/hourly_readings.csv');save_csv(test.reset_index(),p/'data/test_predictions.csv')
    monthly=test.groupby('period').agg(hours=('kwh','size'),actual_kwh=('kwh','sum'),forecast_kwh=('model','sum'),absolute_error=('absolute_error','sum'),baseline_error=('baseline_error','sum'),covered_hours=('covered','sum')).reset_index()
    daily=test.groupby(['day','period']).agg(hours=('kwh','size'),actual_kwh=('kwh','sum'),forecast_kwh=('model','sum'),absolute_error=('absolute_error','sum'),baseline_error=('baseline_error','sum'),covered_hours=('covered','sum')).reset_index()
    expected=frame[frame.index>='2010-01-01'].groupby(frame[frame.index>='2010-01-01'].index.strftime('%Y-%m')).size()
    monthly['calendar_hours']=monthly.period.map(expected);save_csv(monthly,p/'outputs/monthly.csv')
    pattern=test.groupby(['period','hour']).agg(hours=('kwh','size'),mean_kwh=('kwh','mean'),mean_forecast=('model','mean'),mae_kwh=('absolute_error','mean'),covered_hours=('covered','sum')).reset_index();save_csv(pattern,p/'outputs/hour_patterns.csv')
    coverage=hourly.assign(month=hourly.timestamp.dt.strftime('%Y-%m')).groupby('month').agg(calendar_hours=('timestamp','size'),complete_hours=('kwh','count'),measured_minutes=('measured_minutes','sum'),recorded_minutes=('calendar_minutes','sum')).reset_index();save_csv(coverage,p/'outputs/coverage.csv')
    uncertainty=bootstrap_days(test,'baseline','model','kwh','day')
    metrics={'source_rows':raw,'missing_measurements':missing,'hourly_calendar_rows':len(hourly),'complete_hours':int(hourly.kwh.count()),'test_hours':len(test),'test_calendar_hours':int(expected.sum()),'train_rows':len(train),'model_mae_kwh':float(results.iloc[1].mae_kwh),'baseline_mae_kwh':float(results.iloc[0].mae_kwh),'baseline_name':base_name,'improvement_pct':100*(1-results.iloc[1].mae_kwh/results.iloc[0].mae_kwh),'nominal_interval_coverage':.9,'actual_interval_coverage':float(test.covered.mean()),'interval_radius_kwh':radius,'uncertainty':uncertainty}
    js(p/'outputs/metrics.json',metrics)
    js(p/'outputs/evaluation.json',{'split':'2007–2008 training; January–June 2009 configuration and baseline selection; July–December 2009 interval calibration; January–November 2010 untouched testing. Only hours with 60 observed minute readings are scored.','train_rows':len(train),'test_rows':len(test),'selection':'Select lowest validation MAE; freeze fit. Calibrate a nominal 90% symmetric residual band on a separate six-month block. Test dates are never used for selection.','target':'Energy in a complete clock hour, kWh','features':features,'candidates':configs,'forecast_contract':'At midnight, forecast each coming hour using the same hours from prior days; all feature timestamps precede forecast midnight. Earlier test-day data become available to the following forecast day. Assume finalized prior-day measurements available at midnight.','uncertainty':uncertainty,'interval_limitation':'Temporal dependence and distribution shift violate exchangeability; 90% is a nominal calibration target, not a guaranteed future coverage rate.','timezone_limitation':'Source timestamps are naive local clock readings. No UTC or DST disambiguation is supplied; do not infer physical elapsed UTC hours.'})
    quality(p,{'original_count':raw==2075259,'timestamp_unique':not data.timestamp.duplicated().any(),'source_calendar_reconciles':int(hourly.calendar_minutes.sum())==raw,'measured_count_reconciles':int(hourly.measured_minutes.sum())==raw-missing,'complete_hours_only':bool(test.measured_minutes.eq(60).all()),'no_target_imputation':int(eligible.kwh.isna().sum())==0,'split_order':train.index.max()<val.index.min()<calibration.index.min()<test.index.min(),'all_lags_at_least_24_hours':all(int(f.split('_')[-1])>=24 for f in features if f.startswith('lag_')),'finite_predictions':bool(np.isfinite(test.model).all()),'interval_order':bool((test.lower<=test.upper).all()),'monthly_counts_reconcile':int(monthly.hours.sum())==len(test),'energy_conservation':bool(np.isclose(hourly.observed_kwh.sum(),data.Global_active_power.sum()/60))})
    fig,axes=plt.subplots(2,2,figsize=(12,7.8))
    tail=test.iloc[:24*7];axes[0,0].plot(tail.index,tail.kwh,label='Actual',color=INK,lw=1);axes[0,0].plot(tail.index,tail.model,label='Forecast',color=BLUE,lw=1);axes[0,0].fill_between(tail.index,tail.lower,tail.upper,color=BLUE,alpha=.1);axes[0,0].set(title='First observed test week',ylabel='kWh per complete hour');axes[0,0].tick_params(axis='x',rotation=25);axes[0,0].legend(fontsize=8)
    x=np.arange(len(monthly));axes[0,1].plot(x,monthly.covered_hours/monthly.hours*100,'o-',color=BLUE);axes[0,1].axhline(90,color=ORANGE,ls='--');axes[0,1].set_xticks(x,monthly.period.str[5:],rotation=25);axes[0,1].set(title='Interval coverage drifts over time',xlabel='2010 month',ylabel='Observed coverage · %')
    axes[1,0].barh(results.model,results.mae_kwh,color=[TEAL,BLUE]);axes[1,0].set(title='Future-hour forecasting benchmark',xlabel='Mean absolute error · kWh')
    hours=test.groupby('hour')[['kwh','model']].mean();axes[1,1].plot(hours.index,hours.kwh,label='Actual',color=INK);axes[1,1].plot(hours.index,hours.model,label='Forecast',color=BLUE);axes[1,1].set(title='Daily demand profile',xlabel='Local clock hour',ylabel='Mean kWh');axes[1,1].legend(fontsize=8)
    fig.suptitle('Electricity demand planning | One household, 2006–2010',fontsize=17,fontweight='bold');fig.tight_layout();fig.savefig(p/'outputs/chart.png',dpi=150);plt.close(fig)
    js(p/'outputs/dashboard.json',{'kind':'energy','metrics':metrics,'monthly':json.loads(monthly.to_json(orient='records')),'daily':json.loads(daily.to_json(orient='records')),'patterns':json.loads(pattern.to_json(orient='records'))})
    return metrics

if __name__=='__main__': print(json.dumps(run(),indent=2))