#!/bin/sh

# Damon Hart-Davis licenses this file to you
# under the Apache Licence, Version 2.0 (the "Licence");
# you may not use this file except in compliance
# with the Licence. You may obtain a copy of the Licence at
# 
# http://www.apache.org/licenses/LICENSE-2.0
# 
# Unless required by applicable law or agreed to in writing,
# software distributed under the Licence is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the Licence for the
# specific language governing permissions and limitations
# under the Licence.
# 
# Author(s) / Copyright (s): Damon Hart-Davis 2018

# Grid-tied behind-the-meter storage (eg AC-coupled battery) simulator.
# Computes grid flows, and use of storage.
# Can be run with varying levels of storage etc (overriding the defaults)
# to estimate behaviour of each.
#
# Can be run for more than one day's generation data.
#
# Starts with storage empty.
#
# TODO: most of the defaults will be overridable.

# If GRAPHOUT=true then generate output suitable for graphing, eg with gnuplot,
# else generate a textual summary.
#
# Supported values are:
#   * false  no graph output, just a textual summary
#   * daymean  daily mean (per minute) of key parameters
GRAPHOUT=false

# File inputs are simulation of night load and solar input.
# Optionally some scheduled one-off loads such as dishwasher can be provided.

# File for solar PV generation 1-minute (W) samples.
# Can be part day or multiple days, eg a whole month.
# If ending in .gz will be un-GZIPped first.
#
# Sample/format:
#20171201T07:00Z 0
#20171201T08:17Z 1
#20171201T08:18Z 9
#20171201T08:19Z 13
#20171201T08:20Z 14
#20171201T08:21Z 15
#
# Note the date element of the timestamp is ignored.
# The time must be UTC.
#
FPVGENW=data/SunnyBeam/201712.gz

# File for cyclic house loads with 1-minute (W) samples:
#
# Sample/format:
#84
#78
#76
#76
#75
#
FLOADCYCW=img/PV/storesim/fridge-40W-1m-cycle.dat

# Constant/flat household consumption (W).
# 10W typical for 16WW overnight (fans + displays, etc).
FLATW=10

# Storage system self-consumption (W).
# 10W for eg Envoy-S + 1 x Enphase AC Battery inferred.
SSSCW=10

# Effective/usable storage capacity (Wh).
# 1 x Enphase AC Battery ~ 1.1kWh usable.
SCAPWh=1100

# Minimum supported charge/discharge.
# Assumed ~5W for Enphase AC Battery.
SMINCDW=5

# Maximum supported charge/discharge.
# Assumed ~260W for Enphase AC Battery.
SMAXCDW=260

# Storage round-trip efficiency (fraction).
# Full round-trip AC->inverter->battery->inverter->AC.
SRTEFF=0.9

# Previous minute.
# Any missing minutes are filled in, so a run starts at 00:00.
# No fill-in is done if this is an empty string.
# Pad out tail end also to make whole day if not empty.
PREVHHMM=23:59


# Allow overrides of parameters of form PARAM=value.
while [ "$#" -gt 0 ];
do
    ARG=$1
    VALUE="`echo $ARG|awk -F= '{print $NF}'`"
    shift
    case $ARG in
        GRAPHOUT=*) GRAPHOUT=$VALUE;;
        FPVGENW=*) FPVGENW=$VALUE;;
        FLOADCYCW=*) FLOADCYCW=$VALUE;;
        PREVHHMM=*) PREVHHMM=$VALUE;;
        FLATW=*) FLATW=$VALUE;;
        SSSCW=*) SSSCW=$VALUE;;
        SCAPWh=*) SCAPWh=$VALUE;;
        SMINCDW=*) SMINCDW=$VALUE;;
        SMAXCDW=*) SMAXCDW=$VALUE;;
        SRTEFF=*) SRTEFF=$VALUE;;
        *) echo "Unrecognised argument $ARG" 1>&2; exit 1;;
    esac
done


# Check file validity/availability.
if [ ! -s $FPVGENW -o ! -r $FPVGENW ]; then
    echo "PV generation file missing, zero length or not readable: $FPVGENW" 1>&2
    exit 1
fi
if [ ! -s $FLOADCYCW -o ! -r $FLOADCYCW ]; then
    echo "Cyclic house load file missing, zero length or not readable: $FLOADCYCW" 1>&2
    exit 1
fi


# If generating textual summary then dump simulation parameters actually used.
if [ "false" = $GRAPHOUT ]; then
    echo NOT graphing.
    echo "FPVGENW: file input for PV generation: $FPVGENW"
    echo "FLOADCYCW: file input for cyclic house loads: $FLOADCYCW"
    if [ "" != $PREVHHMM ]; then
        echo "PREVHHMM: PV generation padded from after / to: $PREVHHMM"
    fi
    echo "FLATW: constant/flat household consumption (W): $FLATW"
    echo "SSSCW: storage system self-consumption (W): $SSSCW"
    echo "SCAPWh: effective/usable storage capacity (Wh): $SCAPWh"
    echo "SMINCDW: minimum supported charge/discharge (W): $SMINCDW"
    echo "SMAXCDW: maximum supported charge/discharge (W): $SMAXCDW"
    echo "SRTEFF: storage round-trip efficiency (fraction): $SRTEFF"
fi

# Create the pipeline to generate the PV.
# This decompresses and fills in gaps (eg overnight) as needed.
#
# The pipeline output is of the form:
#08:16 0
#08:17 1
#08:17 9
GENGZ=false
GENDECOMP=cat
if [ "gz" = "`echo $FPVGENW|awk -F. '{print $NF}'`" ]; then
    GENGZ=true
    GENDECOMP="gzip -d"
fi
# Input line format:
# Silently skips lines with invalid format.
# (Eg one line in middle of 201712 log is prefixed with NULs.)
#20171201T08:21Z 15
$GENDECOMP < $FPVGENW | awk '
    /^20[0-9]{6}T[0-9]{2}:[0-9]{2}Z [0-9]+$/ {
    HHMM=substr($1, 10, 5);
    # While ++PREVHHMM != HHMM parsed from current record
    # insert filler record with previous watt value.
    while("" != INITPREVHHMM) {
        if("" == PREVHHMM) { PREVHHMM = INITPREVHHMM; }
        split(PREVHHMM, a, ":");
        HH=a[1];
        MM=a[2];
        if(++MM > 59) { MM = 0; if(++HH > 23) { HH = 0; } }
        PREVHHMM=sprintf("%02d:%02d", HH, MM);
        if(PREVHHMM == HHMM) { break; }
        print PREVHHMM, PREVW;
        }
    # Then print current record and update PREVHHMM to it.
    print HHMM, $2
    PREVHHMM=HHMM;
    PREVW=$2
    }
    END {
    # If start was potentally padded,
    # then pad at end to make full day (generally with zeros).
    while("" != INITPREVHHMM) {
        split(PREVHHMM, a, ":");
        HH=a[1];
        MM=a[2];
        if(++MM > 59) { MM = 0; if(++HH > 23) { HH = 0; } }
        PREVHHMM=sprintf("%02d:%02d", HH, MM);
        print PREVHHMM, PREVW;
        if(PREVHHMM == INITPREVHHMM) { break; }
        }
    }
    ' INITPREVHHMM=$PREVHHMM PREVW=0 | \
awk '{
    # Process each record in turn.
    # Count records
    ++mins;
    HHMM=$1;
    ++count[HHMM];
    # 0/+ve generation in watts; TgenWm is total in watt-minutes.
    genW = $2;
    TgenWm += genW;
    # 0/+ve load in watts; TloadWm is total in watt-minutes.
    loadW = FLATW + SSSCW;
    # Add in cyclic load.
    if("" != FLOADCYCW) {
        if((getline clW < FLOADCYCW) <= 0) {
            # If cyclic load file has finished, close and reopen.
            close(FLOADCYCW);
            getline clW < FLOADCYCW;
            }
        else {
#print "clW", clW;
            }
        if(clW > 0) { loadW += clW; }
        }
    TloadWm += loadW;
    # Raw net flow (without batteries) in watts OUT to grid.
    netflowW = genW - loadW;
    TnetflowWm += netflowW;
    # Absoute raw net flow (without batteries) in watts.
    absnetflowW = netflowW;
    if(absnetflowW < 0) { absnetflowW = -absnetflowW; }
    TabsnetflowWm += absnetflowW;
    # Total raw net imports (without battery) from grid (-ve) (Wm).
    if(netflowW < 0) { TrawimportsWm += netflowW; }
    # Total raw net exports (without battery) to grid (+ve) (Wm).
    if(netflowW > 0) { TrawexportsWm += netflowW; }
    # Adjustment charge flow into battery (+ve) or discharge (after losses).
    adjgridflowWm = 0;
    # Conservatively, all losses are taken as happening at discharge.
    # Battery running stored charge is in battChargeWm (initally 0).
    # Battery charge?
    # Flow into (and thus through) battery (Wm).
    battthroughputW = 0;
    # Generation - load has to be +ve and at least SMINCDW.
    # There also has to be space in the battery.
    if(netflowW > SMINCDW) {
        battSpaceWm = (60*SCAPWh) - battChargeWm;
#print "battSpaceWm", battSpaceWm;
        if(battSpaceWm > 0) {
            # Only flow above the low threshold and up to max is corrected.
            correctableFlowW = netflowW - SMINCDW;
            if(correctableFlowW > SMAXCDW) { correctableFlowW = SMAXCDW; }
            # Treat all losses as happening at discharge.
            # This treates battery capacity conservatively.
            if(battSpaceWm > correctableFlowW) {
                f = correctableFlowW;
                } else {
                f = battSpaceWm;
                }
            adjgridflowWm = f;
            battChargeWm += f;
            battthroughputW = f;
            }
        }
    # Battery discharge?
    # Generation - load has to be -ve and at/below -SMINCDW.
    # There also has to be charge in the battery.
    else if(netflowW < -SMINCDW) {
        if(battChargeWm > 0) {
            # Only flow above the low threshold and up to max is corrected.
            correctableFlowW = -netflowW - SMINCDW;
            if(correctableFlowW > SMAXCDW) { correctableFlowW = SMAXCDW; }
            # Treat all losses as happening at discharge.
            # This treates battery capacity conservatively.
            availableBattChargeWm = SRTEFF * battChargeWm;
            if(availableBattChargeWm < correctableFlowW) {
                f = availableBattChargeWm;
                } else {
                f = correctableFlowW;
                }
            adjgridflowWm = -f;
            battChargeWm -= (f / SRTEFF);
            }
        }
    chargeWm[HHMM] += battChargeWm;
#print "battChargeWm", battChargeWm, "adjgridflowWm", adjgridflowWm, HHMM;
    # Sum battery charge (flow).
    TbattthroughputWm += battthroughputW;
    # Sum absolute avoided grid flows (Wm).
#print "adjgridflowWm", adjgridflowWm;
    if(adjgridflowWm > 0) { TadjgridflowWm += adjgridflowWm;}
    else { TadjgridflowWm -= adjgridflowWm; }
    # Net final flow (with batteries) OUT to grid in watts.
    finalgridflowW = genW - loadW - adjgridflowWm;
#print "finalgridflowW", finalgridflowW;
    flow[HHMM] += finalgridflowW;
    # Total raw net imports (with battery) from grid (-ve) (Wm).
    if(finalgridflowW < 0) { TfinalimportsWm += finalgridflowW; }
    # Total raw net exports (with battery) to grid (+ve) (Wm).
    if(finalgridflowW > 0) { TfinalexportsWm += finalgridflowW; }
    }
    END {
        # Print totals (if not creating a graph).
        if("false" == GRAPHOUT) {
            print "*** RESULTS";
            print "Sim mins", mins, \
                  "| (Wm)", \
                  "TgenWm", TgenWm, \
                  "TloadWm", TloadWm, \
                  "TnetflowWm", TnetflowWm, "TabsnetflowWm", TabsnetflowWm, \
                  "TbattthroughputWm", TbattthroughputWm, "TadjgridflowWm", TadjgridflowWm, \
                  "TrawimportsWm", TrawimportsWm, \
                  "TrawexportsWm", TrawexportsWm, \
                  "TfinalimportsWm", TfinalimportsWm, \
                  "TfinalexportsWm", TfinalexportsWm;
            days = mins / 1440;
            cf = 1 / (60000 * days); # Conversion factor to kWh/d.
            print "Sim days", days, \
                  "| (kWh/d)", \
                  "Tgen", TgenWm * cf, \
                  "Tload", TloadWm * cf, \
                  "Tnetflow", TnetflowWm*cf, "Tabsnetflow", TabsnetflowWm*cf, \
                  "Tbattthroughput", TbattthroughputWm*cf, "Tadjgridflow", TadjgridflowWm*cf,\
                  "Trawimports", TrawimportsWm*cf, \
                  "Trawexports", TrawexportsWm*cf, \
                  "Tfinalimports", TfinalimportsWm*cf, \
                  "Tfinalexports", TfinalexportsWm*cf;
            }
        else if("daymean" == GRAPHOUT) {
            # "Daily mean" by minute graph.
            # Designed to be consumed by gnuplot or similar.
            # Each line starts "G daymean ".
            for(HH = 0; HH < 24; ++HH) {
                for(MM = 0; MM < 60; ++MM) {
                    key = sprintf("%02d:%02d", HH, MM);
                    n = count[key];
                    f = flow[key] / n;
                    c = chargeWm[key] / ((60*SCAPWh)/100 * n);
                    printf("G daymean HH:MM %s n %d flow|W %d SoC|%% %d\n",
                        key, n, f, c);
                    }
                }
            }
    }' \
    GRAPHOUT=$GRAPHOUT \
    FLATW=$FLATW SSSCW=$SSSCW SMINCDW=$SMINCDW SMAXCDW=$SMAXCDW SCAPWh=$SCAPWh \
    SRTEFF=$SRTEFF \
    FLOADCYCW=$FLOADCYCW