feat(application): integrate measured-load ingestion training and planner source
This commit is contained in:
@@ -0,0 +1,199 @@
|
||||
from __future__ import annotations
|
||||
from math import sqrt
|
||||
from uuid import uuid4
|
||||
import numpy as np
|
||||
from scipy.optimize import Bounds, LinearConstraint, milp
|
||||
from scipy.sparse import coo_matrix
|
||||
from .domain import Battery, Limits, Step, month_key, number, quarter_start, utc
|
||||
from .peak_policy import basis_record, RestMonthOutlook
|
||||
|
||||
class Model:
|
||||
def __init__(self):
|
||||
self.lower,self.upper,self.cost,self.integer=[],[],[],[]
|
||||
self.rows,self.row_lo,self.row_hi=[],[],[]
|
||||
def variable(self, lower=0., upper=np.inf, cost=0., integer=False):
|
||||
index=len(self.lower)
|
||||
self.lower.append(lower);self.upper.append(upper);self.cost.append(cost);self.integer.append(int(integer))
|
||||
return index
|
||||
def constraint(self, coefficients, lower=-np.inf, upper=np.inf):
|
||||
self.rows.append(coefficients);self.row_lo.append(lower);self.row_hi.append(upper)
|
||||
def matrices(self):
|
||||
rr,cc,vv=[],[],[]
|
||||
for row,coefficients in enumerate(self.rows):
|
||||
for col,value in coefficients.items():
|
||||
rr.append(row);cc.append(col);vv.append(value)
|
||||
matrix=coo_matrix((vv,(rr,cc)),shape=(len(self.rows),len(self.lower))).tocsr()
|
||||
return matrix,np.array(self.row_lo),np.array(self.row_hi)
|
||||
|
||||
from .battery_model import BatteryModel
|
||||
|
||||
class AssetRegistry:
|
||||
"""Reviewed asset classes only; future EV/thermal models add their own constraints."""
|
||||
def __init__(self):self.models={Battery:BatteryModel}
|
||||
def register(self,asset_type,implementation):
|
||||
if asset_type in self.models:raise ValueError('Asset model already registered')
|
||||
self.models[asset_type]=implementation
|
||||
def build(self,model,asset,steps,direction):
|
||||
if type(asset) not in self.models:raise ValueError('Unsupported asset model')
|
||||
return self.models[type(asset)].build(model,asset,steps,direction)
|
||||
|
||||
def _error(reason,status='invalid_inputs'):
|
||||
return {'schemaVersion':2,'status':status,'executable':False,'points':[],'reason':str(reason)}
|
||||
|
||||
def _validate(steps,assets,at,past,observed_peaks,peak_prices,limits):
|
||||
if not steps or len(steps)>576:raise ValueError('Need 1..576 steps')
|
||||
for i,step in enumerate(steps):
|
||||
start=utc(step.start)
|
||||
if type(step.seconds) is not int or not 1<=step.seconds<=300:raise ValueError('Interval duration must be 1..300 seconds')
|
||||
if start.microsecond or step.end.second or step.end.microsecond or step.end.minute%5:raise ValueError('Intervals must end on a 5-minute boundary')
|
||||
if i and (step.seconds!=300 or start.second or start.minute%5):raise ValueError('Only first interval may be shortened')
|
||||
if i and start!=steps[i-1].end:raise ValueError('Missing/duplicated/overlapping interval')
|
||||
number(step.base_load_w,'base load',0);number(step.pv_w,'PV',0);number(step.external_w,'external flow')
|
||||
if not step.import_price or not step.export_price or not step.import_price.known_at(at) or not step.export_price.known_at(at):raise ValueError('Unknown or unpublished prices')
|
||||
limits.import_limit(start)
|
||||
if utc(at)!=utc(steps[0].start):raise ValueError('First interval must start at decision time')
|
||||
if steps[-1].end!=quarter_start(steps[-1].end):raise ValueError('End horizon on complete billing quarter')
|
||||
if limits.export_w is not None:number(limits.export_w,'export limit',0)
|
||||
ids=set()
|
||||
for asset in assets:
|
||||
asset.validate(at)
|
||||
if asset.asset_id in ids:raise ValueError('Duplicate asset_id')
|
||||
ids.add(asset.asset_id)
|
||||
q=quarter_start(steps[0].start);elapsed=int((utc(steps[0].start)-q).total_seconds())
|
||||
if any(key!=q for key in past):raise ValueError('Only elapsed energy of first quarter permitted')
|
||||
if not elapsed and q in past and (past[q].import_kwh!=0 or past[q].measured_seconds!=0):raise ValueError('No elapsed energy at quarter boundary')
|
||||
if elapsed:
|
||||
if q not in past or past[q].measured_seconds!=elapsed:raise ValueError('Actual elapsed quarter import energy missing')
|
||||
number(past[q].import_kwh,'quarter import energy',0)
|
||||
for month in {month_key(s.start) for s in steps}:
|
||||
if month not in observed_peaks or month not in peak_prices:raise ValueError(f'Measured peak state or tariff missing for {month}')
|
||||
number(observed_peaks[month],'measured peak',0);number(peak_prices[month],'peak tariff',0)
|
||||
|
||||
def optimize(steps,batteries,*,at,limits=None,observed_peaks=None,peak_prices=None,quarter_history=None,config_revision=1,family='3',timeout_seconds=30.,asset_registry=None,peak_context=None,peak_outlooks=None):
|
||||
"""Pure MILP. Peak state is measured, never a configured cap. No device calls."""
|
||||
limits=limits or Limits();observed_peaks=observed_peaks or {};peak_prices=peak_prices or {}
|
||||
past={utc(k):v for k,v in (quarter_history or {}).items()}
|
||||
try:
|
||||
_validate(steps,batteries,at,past,observed_peaks,peak_prices,limits)
|
||||
number(timeout_seconds,'solver timeout',.01,600)
|
||||
contexts={}
|
||||
for month in {month_key(s.start) for s in steps}:
|
||||
if peak_context is None:
|
||||
contexts[month]={'kw':observed_peaks[month],'quality':'verified','source':'legacy_verified_contract','observedAt':utc(at).isoformat()}
|
||||
else:
|
||||
contexts[month]=basis_record(peak_context[month],month,at,allow_estimates=True)
|
||||
if abs(contexts[month]['kw']-observed_peaks[month])>1e-8:raise ValueError('Peak context differs from numerical basis')
|
||||
outlooks=peak_outlooks or {}
|
||||
if set(outlooks)-set(contexts):raise ValueError('Outlook for month outside horizon')
|
||||
for month,outlook in outlooks.items():
|
||||
if not isinstance(outlook,RestMonthOutlook) or outlook.month!=month:raise ValueError('Invalid peak outlook')
|
||||
outlook.validate(at,steps[-1].end)
|
||||
except (ValueError,TypeError,KeyError,AttributeError) as exc:return _error(exc)
|
||||
model=Model();direction=[model.variable(0,1,integer=True) for _ in steps];registry=asset_registry or AssetRegistry()
|
||||
try:handles={b.asset_id:registry.build(model,b,steps,direction) for b in batteries}
|
||||
except ValueError as exc:return _error(exc)
|
||||
total_charge=sum(b.max_charge_w for b in batteries)/1000
|
||||
total_discharge=sum(b.max_discharge_w for b in batteries)/1000
|
||||
imp,exp,curtail=[],[],[];quarters={}
|
||||
for i,step in enumerate(steps):
|
||||
residual=step.residual_w/1000;dt=step.seconds/3600
|
||||
imax=max(0.,(step.base_load_w+step.external_w)/1000)+total_charge
|
||||
emax=max(0.,-residual)+total_discharge;limit=limits.import_limit(step.start)
|
||||
if limit is not None:imax=min(imax,limit/1000)
|
||||
if limits.export_w is not None:emax=min(emax,limits.export_w/1000)
|
||||
pi=model.variable(0,imax,step.import_price.chf_kwh*dt)
|
||||
pe=model.variable(0,emax,-step.export_price.chf_kwh*dt)
|
||||
pc=model.variable(0,step.pv_w/1000)
|
||||
imp.append(pi);exp.append(pe);curtail.append(pc)
|
||||
gm=model.variable(0,1,integer=True)
|
||||
model.constraint({pi:1,gm:-imax},upper=0);model.constraint({pe:1,gm:emax},upper=emax)
|
||||
balance={pi:1,pe:-1,pc:-1};pv_only={}
|
||||
for b in batteries:
|
||||
h=handles[b.asset_id];balance[h['charge'][i]]=-1;balance[h['discharge'][i]]=1
|
||||
if not b.grid_charging:pv_only[h['charge'][i]]=1
|
||||
model.constraint(balance,residual,residual)
|
||||
if pv_only:
|
||||
model.constraint(pv_only,upper=max(0.,-residual))
|
||||
no_grid_max=sum(b.max_charge_w for b in batteries if not b.grid_charging)/1000
|
||||
model.constraint({**pv_only,gm:no_grid_max},upper=no_grid_max)
|
||||
quarters.setdefault(quarter_start(step.start),[]).append(i)
|
||||
months=sorted({month_key(s.start) for s in steps})
|
||||
peaks={m:model.variable(observed_peaks[m]) for m in months}
|
||||
for m in months:
|
||||
if m in outlooks:outlooks[m].add_to_model(model,peaks[m],observed_peaks[m],peak_prices[m])
|
||||
else:model.cost[peaks[m]]=peak_prices[m]
|
||||
for quarter,positions in quarters.items():
|
||||
coefficients={imp[i]:steps[i].seconds/900 for i in positions};coefficients[peaks[month_key(quarter)]]=-1
|
||||
used=past[quarter].import_kwh if quarter in past else 0.
|
||||
model.constraint(coefficients,upper=-used/.25)
|
||||
matrix,row_lo,row_hi=model.matrices()
|
||||
try:
|
||||
result=milp(np.array(model.cost),integrality=np.array(model.integer),bounds=Bounds(model.lower,model.upper),constraints=LinearConstraint(matrix,row_lo,row_hi),options={'time_limit':float(timeout_seconds),'mip_rel_gap':1e-4})
|
||||
except Exception as exc:return _error(f'Solver exception: {type(exc).__name__}','solver_error')
|
||||
x=result.x
|
||||
if result.status not in (0,1) or x is None:return _error(result.message,'no_feasible_plan')
|
||||
if not np.isfinite(x).all():return _error('Non-finite solver result','validation_failed')
|
||||
ax=matrix@x;tol=1e-6
|
||||
if (np.any(x<np.array(model.lower)-tol) or np.any(x>np.array(model.upper)+tol) or np.any(ax<row_lo-tol) or np.any(ax>row_hi+tol) or any(abs(x[i]-round(x[i]))>tol for i,flag in enumerate(model.integer) if flag)):
|
||||
return _error('Solver incumbent violates constraints','validation_failed')
|
||||
running_peaks=dict(observed_peaks);baseline_peaks=dict(observed_peaks);peak_deltas={};baseline_peak_deltas={}
|
||||
for quarter,positions in quarters.items():
|
||||
used=past[quarter].import_kwh if quarter in past else 0.;m=month_key(quarter)
|
||||
demand=(used+sum(x[imp[i]]*steps[i].seconds/3600 for i in positions))/.25
|
||||
baseline=(used+sum(max(0.,steps[i].residual_w)*steps[i].seconds/3600000 for i in positions))/.25
|
||||
before=running_peaks[m];running_peaks[m]=max(before,demand);peak_deltas[positions[-1]]=(running_peaks[m]-before)*peak_prices[m]
|
||||
before=baseline_peaks[m];baseline_peaks[m]=max(before,baseline);baseline_peak_deltas[positions[-1]]=(baseline_peaks[m]-before)*peak_prices[m]
|
||||
points=[];cumulative_energy=cumulative_cash=cumulative_baseline_energy=cumulative_baseline_cash=cumulative_throughput=0.
|
||||
for i,step in enumerate(steps):
|
||||
dt=step.seconds/3600;p_import=max(0.,float(x[imp[i]]));p_export=max(0.,float(x[exp[i]]))
|
||||
energy_cost=(p_import*step.import_price.chf_kwh-p_export*step.export_price.chf_kwh)*dt
|
||||
baseline_import=max(0.,step.residual_w)/1000;baseline_export=max(0.,-step.residual_w)/1000
|
||||
if limits.export_w is not None:baseline_export=min(baseline_export,limits.export_w/1000)
|
||||
baseline_cost=(baseline_import*step.import_price.chf_kwh-baseline_export*step.export_price.chf_kwh)*dt
|
||||
targets,soc_end={},{};throughput_cost=0.
|
||||
for b in batteries:
|
||||
h=handles[b.asset_id];targets[b.asset_id]=float((x[h['charge'][i]]-x[h['discharge'][i]])*1000)
|
||||
soc_end[b.asset_id]=float(100*x[h['energy'][i+1]]/b.capacity_kwh)
|
||||
throughput_cost+=float((x[h['charge'][i]]+x[h['discharge'][i]])*dt*b.throughput_chf_kwh)
|
||||
cumulative_energy+=energy_cost;cumulative_cash+=energy_cost+peak_deltas.get(i,0.)
|
||||
cumulative_baseline_energy+=baseline_cost;cumulative_baseline_cash+=baseline_cost+baseline_peak_deltas.get(i,0.)
|
||||
cumulative_throughput+=throughput_cost;battery_w=sum(targets.values())
|
||||
if battery_w>1:intent='gridCharge' if p_import>.001 else 'pvCharge'
|
||||
elif battery_w< -1:intent='export' if p_export>.001 else 'discharge'
|
||||
else:intent='hold'
|
||||
points.append({'time':utc(step.start).isoformat(),'validUntil':step.end.isoformat(),
|
||||
'baselineGridW':float(step.residual_w),'gridTargetW':(p_import-p_export)*1000,
|
||||
'batteryTargetW':battery_w,'assetTargetsW':targets,'socEndPercent':soc_end,
|
||||
'pvCurtailmentW':float(x[curtail[i]]*1000),'intent':intent,
|
||||
'importLimitW':limits.import_limit(step.start),'exportLimitW':limits.export_w,
|
||||
'importPriceChfKwh':step.import_price.chf_kwh,'exportPriceChfKwh':step.export_price.chf_kwh,
|
||||
'energyCostChf':energy_cost,'additionalPeakCostChf':peak_deltas.get(i,0.),
|
||||
'throughputCostChf':throughput_cost,'cumulativeEnergyCostChf':cumulative_energy,
|
||||
'cumulativeCashCostChf':cumulative_cash,'cumulativeBaselineEnergyCostChf':cumulative_baseline_energy,
|
||||
'cumulativeBaselineCashCostChf':cumulative_baseline_cash,'baselineAdditionalPeakCostChf':baseline_peak_deltas.get(i,0.)})
|
||||
end_value=sum(float(x[handles[b.asset_id]['energy'][-1]])*b.terminal_value_chf_kwh for b in batteries)
|
||||
def planning_cost(chosen):
|
||||
return sum(outlooks[m].incremental_cost(observed_peaks[m],chosen[m],peak_prices[m]) if m in outlooks
|
||||
else max(0.,chosen[m]-observed_peaks[m])*peak_prices[m] for m in months)
|
||||
planning_peak=planning_cost(running_peaks)
|
||||
full_peak=cumulative_cash-cumulative_energy
|
||||
estimated=any(c['quality']=='estimated' for c in contexts.values())
|
||||
return {'schemaVersion':2,'planId':str(uuid4()),'configRevision':config_revision,'sourceFamily':family,
|
||||
'status':'optimal' if result.status==0 else 'feasible_time_limit','executable':True,
|
||||
'generatedAt':utc(at).isoformat(),'validFrom':points[0]['time'],'validUntil':points[-1]['validUntil'],
|
||||
'intervalMinutes':5,'points':points,'solverGap':float(result.mip_gap) if getattr(result,'mip_gap',None) is not None else None,
|
||||
'measuredPeaksKw':{m:observed_peaks[m] for m in months if contexts[m]['quality']=='verified'},
|
||||
'peakBasisKw':{m:observed_peaks[m] for m in months},'peakBasis':contexts,
|
||||
'peakCostIsEstimate':estimated,'planningPeakCostChf':planning_peak,
|
||||
'restMonthAdjustmentChf':planning_peak-full_peak,
|
||||
'peakScenarioOutlook':{m:o.as_dict() for m,o in outlooks.items()},
|
||||
'plannedPeaksKw':{m:running_peaks[m] for m in months},
|
||||
'additionalPeakCostChf':cumulative_cash-cumulative_energy,'energyCostChf':cumulative_energy,'cashCostChf':cumulative_cash,
|
||||
'baselineEnergyCostChf':cumulative_baseline_energy,'baselineCashCostChf':cumulative_baseline_cash,
|
||||
'baselineAdditionalPeakCostChf':sum(baseline_peak_deltas.values()),'baselinePeaksKw':baseline_peaks,
|
||||
'throughputCostChf':cumulative_throughput,'terminalValueChf':end_value,
|
||||
'objectiveChf':cumulative_energy+planning_peak+cumulative_throughput-end_value,
|
||||
'baselinePlanningPeakCostChf':planning_cost(baseline_peaks),
|
||||
'cashCostMeaning':'Projected horizon energy plus full incremental monthly tariff, relative to the labelled peak basis; not an invoice',
|
||||
'terminalMinSocPercent':{b.asset_id:handles[b.asset_id]['terminal_min_percent'] for b in batteries},
|
||||
'peakOutlook':'rest_month_scenarios' if outlooks else 'full_incremental_tariff'}
|
||||
Reference in New Issue
Block a user