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(xnp.array(model.upper)+tol) or np.any(axrow_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'}