V2G Microgrid

import numpy as np
from scipy.optimize import linprog

T = 24                                         # 1-hour steps, day-ahead
load = np.array([60,55,52,50,50,55,70,95,120,135,140,142,140,138,135,130,120,110,105,100,90,80,70,65.])
pv   = 1.5*np.array([0,0,0,0,0,2,15,40,70,95,112,120,118,105,85,55,25,6,0,0,0,0,0,0.])
tou  = np.array([4]*6 + [6]*4 + [7]*8 + [10]*4 + [6]*2, float)   # INR/kWh import
feed = 3.0                                     # INR/kWh export
G, eta = 150., .95                             # grid limit (kW), one-way efficiency
BC, BP = 200., 50.                             # BESS kWh, kW
NEV, EC, EP = 10, 40., 7.                      # EV count, kWh, kW each
ARR, DEP, S0, SREQ = 10, 20, .40, .80          # fleet plugged 10:00-20:00
EE, EPM = NEV*EC, NEV*EP
wear_b, wear_e = 1.5, 2.0                      # INR per kWh throughput
DC = 40.                                       # demand charge, INR per kW of daily peak import

# variables: grid-in, grid-out, BESS ch/dis, EV ch/dis, PV curtail  (7*T) + SOC_b, SOC_ev (T+1 each)
n = 7*T + 2*(T+1) + 1; v = lambda k,t: k*T + t
sb = lambda t: 7*T + t;  se = lambda t: 8*T + 1 + t
PK = n - 1
c = np.zeros(n); c[PK] = DC
for t in range(T):
    c[v(0,t)], c[v(1,t)] = tou[t], -feed
    c[v(2,t)] = c[v(3,t)] = wear_b/2
    c[v(4,t)] = c[v(5,t)] = wear_e/2
A, b = [], []
for t in range(T):
    r = np.zeros(n); r[[v(0,t), v(3,t), v(5,t)]] = 1; r[[v(1,t), v(2,t), v(4,t), v(6,t)]] = -1
    A.append(r); b.append(load[t] - pv[t])                             # power balance
    r = np.zeros(n); r[sb(t+1)], r[sb(t)] = 1, -1; r[v(2,t)], r[v(3,t)] = -eta/BC, 1/(eta*BC)
    A.append(r); b.append(0)                                           # BESS energy
    r = np.zeros(n); r[se(t+1)], r[se(t)] = 1, -1; r[v(4,t)], r[v(5,t)] = -eta/EE, 1/(eta*EE)
    A.append(r); b.append(0)                                           # EV energy
Aub = np.zeros((T,n)); bub = np.zeros(T)
for t in range(T): Aub[t,v(0,t)], Aub[t,PK] = 1, -1              # grid import <= peak variable
bd = [(0,0)]*n; bd[PK] = (0,G)
for t in range(T):
    plug = ARR <= t < DEP
    bd[v(0,t)] = bd[v(1,t)] = (0,G)
    bd[v(2,t)] = bd[v(3,t)] = (0,BP)
    bd[v(4,t)] = bd[v(5,t)] = (0, EPM if plug else 0)
    bd[v(6,t)] = (0, pv[t])
for t in range(T+1):
    bd[sb(t)] = (.2,.9)
    bd[se(t)] = (.2,.95) if ARR <= t <= DEP else (0,1)
bd[sb(0)] = (.5,.5); bd[se(ARR)] = (S0,S0); bd[se(DEP)] = (SREQ,.95)
A.append(np.eye(n)[sb(T)] - np.eye(n)[sb(0)]); b.append(0)             # BESS ends where it started

r = linprog(c, A_ub=Aub, b_ub=bub, A_eq=np.array(A), b_eq=b, bounds=bd, method="highs"); x = r.x
gi, ge, bc, bdis, ec, ed = (x[k*T:(k+1)*T] for k in range(6))
cost = np.sum(tou*gi - feed*ge) + DC*x[PK]
wear = wear_b/2*np.sum(bc+bdis) + wear_e/2*np.sum(ec+ed)

# baseline: no BESS dispatch, EVs charge flat-out on arrival
need, ev0 = (SREQ-S0)*EE/eta, np.zeros(T)
for t in range(ARR, DEP): ev0[t] = min(EPM, need); need -= ev0[t]
net = load + ev0 - pv
base = np.sum(tou*np.maximum(net,0) - feed*np.maximum(-net,0)) + DC*np.maximum(net,0).max()
print(f"Baseline cost        : INR {base:8.1f}   peak import {np.maximum(net,0).max():6.1f} kW")
print(f"Optimised (energy+DC): INR {cost:8.1f}   peak import {gi.max():6.1f} kW")
print(f"Battery wear cost    : INR {wear:8.1f}")
print(f"Net saving           : INR {base-cost-wear:8.1f}  ({100*(base-cost-wear)/base:.1f} %)")
print(f"EV energy exported   : {ed.sum():.0f} kWh   BESS discharged: {bdis.sum():.0f} kWh   PV curtailed: {x[6*T:7*T].sum():.0f} kWh")

Baseline cost        : INR  12701.8   peak import  105.0 kW
Optimised (energy+DC): INR  10749.8   peak import   70.1 kW
Battery wear cost    : INR    416.2
Net saving           : INR   1535.7  (12.1 %)
EV energy exported   : 57 kWh   BESS discharged: 81 kWh   PV curtailed: 0 kWh