simulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from matplotlib import pyplot\n",
    "import json\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def load_route(route, reverse=False):\n",
    "    if reverse:\n",
    "        data_file = 'route_data/{}_data.json'.format(route)\n",
    "    else:\n",
    "        data_file = 'route_data/{}_reverse_data.json'.format(route)\n",
    "    return json.load(open(data_file, 'r'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define simulation constants\n",
    "\n",
    "ROUTES = [10, 11, 15, 17, 81, 82]\n",
    "\n",
    "# physical constants\n",
    "GRAVITY              = 9.81  # [m/s^2]\n",
    "FUEL_DENSITY         = 35800 # [kJ/L]\n",
    "AIR_DENSITY          = 1.225 # [kg/m^3]\n",
    "\n",
    "# bus specifications\n",
    "MASS                 = 9500  # [kg]\n",
    "CROSS_AREA           = 7.63  # [m^2]\n",
    "\n",
    "# model parameters\n",
    "parameters = {\n",
    "    'DRAG_COEFFICIENT'     : 0.70,  # []\n",
    "    'ROLLING_COEFFICIENT'  : 0.012, # []\n",
    "    'ENGINE_EFFICIENCY'    : 0.35,  # []\n",
    "    'BRAKING_EFFICIENCY'   : 0.50,  # []\n",
    "    'CHARGING_EFFICIENCY'  : 0.65,  # []\n",
    "    'CAPACITY'             : 44600, # [kj]\n",
    "    'MAX_CHARGE_RATE'      : 37,    # [kW]\n",
    "    'MAX_DISCHARGE_RATE'   : 99,    # [kW]\n",
    "    'ENGINE_CHARGE_RATE'   : 1,     # [kW]\n",
    "    'BRAKING_ACCEL'        : 1.5,   # [m/s^2]\n",
    "    'FORWARD_ACCEL'        : 1.1,   # [m/s^2]\n",
    "    'IDLE_FUEL_USE'        : 3.93,  # [L/h]\n",
    "    'STOP_TIME'            : 12,    # [s]\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# run simulation\n",
    "\n",
    "def run_simulation(route, params, hybrid, initial_battery, time_limit, dt, reverse, winter):\n",
    "    path = load_route(route, reverse)\n",
    "    \n",
    "    position = 0\n",
    "    velocity = 0\n",
    "    time_stopped = 0\n",
    "    fuel_consumed = 0\n",
    "    battery = initial_battery\n",
    "    total_regen = 0\n",
    "\n",
    "    output = []\n",
    "    stop_times = []\n",
    "    last_stop = None\n",
    "\n",
    "    for t in np.arange(0, time_limit * 60, dt):\n",
    "        if position > path[-1]['dist']:\n",
    "            break\n",
    "            \n",
    "        index = len([p for p in path if p['dist'] < position])\n",
    "        next_point = path[index]\n",
    "        last_point = path[index - 1]\n",
    "        \n",
    "        stop_offset = 3\n",
    "        if next_point['stop'] and next_point['dist'] - position < stop_offset*1.1 and last_stop != index - 1:\n",
    "            stop_times.append(t / 60)\n",
    "            last_stop = index - 1\n",
    "\n",
    "        speed_limit = 15.6\n",
    "        gradient = last_point['gradient']\n",
    "        elevation = (last_point['elevation'] if index == 0\n",
    "                else last_point['elevation'] + (position - last_point['dist']) * gradient)\n",
    "\n",
    "        future_stops = [p['dist'] for p in path[index:] if p['stop']]\n",
    "        dist_to_stop = 1e9 if len(future_stops) == 0 else future_stops[0] - position\n",
    "        stopping_time = velocity / params['BRAKING_ACCEL']\n",
    "        stopping_dist = 0.5 * velocity * stopping_time\n",
    "        \n",
    "        # stop and wait if at a stop\n",
    "        if dist_to_stop < stop_offset and time_stopped < params['STOP_TIME']:\n",
    "            if velocity > params['BRAKING_ACCEL'] * dt:\n",
    "                acceleration = -params['BRAKING_ACCEL']\n",
    "            else:\n",
    "                acceleration = -velocity * dt\n",
    "            time_stopped += dt\n",
    "\n",
    "        # slow down if approaching a stop\n",
    "        elif stopping_dist + stop_offset >= dist_to_stop and time_stopped < params['STOP_TIME']:\n",
    "            acceleration = -params['BRAKING_ACCEL']\n",
    "\n",
    "        # speed up if not approaching a stop\n",
    "        elif velocity < speed_limit:\n",
    "            acceleration = params['FORWARD_ACCEL']\n",
    "\n",
    "        # maintain speed\n",
    "        else:\n",
    "            acceleration = 0\n",
    "            \n",
    "        if dist_to_stop > stop_offset:\n",
    "            time_stopped = 0\n",
    "\n",
    "        if winter and gradient > 0.05:\n",
    "            slippage_factor = 1 / 0.8\n",
    "        else:\n",
    "            slippage_factor = 1\n",
    "            \n",
    "        power = 0.001 * (\n",
    "                MASS * GRAVITY * gradient * velocity + # gravitational\n",
    "                MASS * velocity * acceleration * slippage_factor + # accelerational\n",
    "                0.5 * params['DRAG_COEFFICIENT'] * AIR_DENSITY * CROSS_AREA * velocity ** 3 + # air drag\n",
    "                MASS * GRAVITY * params['ROLLING_COEFFICIENT'] * velocity # rolling drag\n",
    "        )\n",
    "\n",
    "        if hybrid:\n",
    "            real_drate = min(battery / dt, params['MAX_DISCHARGE_RATE'])\n",
    "            real_crate = min((params['CAPACITY'] - battery) / dt, params['MAX_CHARGE_RATE'])\n",
    "            battery_power = min(real_drate, max(-real_crate, power))\n",
    "\n",
    "            if battery >= params['CAPACITY']:\n",
    "                battery_power = max(battery_power, 0)\n",
    "\n",
    "            if battery <= 0:\n",
    "                battery_power = min(battery_power, 0)\n",
    "\n",
    "            if battery_power > 0: # discharging\n",
    "                battery -= battery_power * dt\n",
    "            else: # charging\n",
    "                battery -= battery_power * params['BRAKING_EFFICIENCY'] * dt\n",
    "                total_regen -= battery_power * params['BRAKING_EFFICIENCY'] * dt\n",
    "\n",
    "            if power > battery_power:\n",
    "                engine_power = power - battery_power\n",
    "            else:\n",
    "                engine_power = 0\n",
    "\n",
    "            if engine_power > 0 and battery < params['CAPACITY']:\n",
    "                engine_power += params['ENGINE_CHARGE_RATE']\n",
    "                battery += params['ENGINE_CHARGE_RATE'] * params['CHARGING_EFFICIENCY'] * dt\n",
    "\n",
    "        else:\n",
    "            battery_power = 0\n",
    "            engine_power = max(power, 0)\n",
    "            \n",
    "        idle_drain = params['IDLE_FUEL_USE'] / 3600 * FUEL_DENSITY * params['ENGINE_EFFICIENCY']\n",
    "        engine_power = max(idle_drain, engine_power)\n",
    "        \n",
    "        fuel_consumed += engine_power / params['ENGINE_EFFICIENCY'] / FUEL_DENSITY * dt\n",
    "\n",
    "        velocity += acceleration * dt\n",
    "        position += velocity * dt\n",
    "\n",
    "        output.append([\n",
    "                ['Time [min]', t / 60], \n",
    "                ['Velocity [m/s]', velocity], \n",
    "                ['Elevation [m]', elevation],\n",
    "                ['Gradient', gradient], \n",
    "                ['Power [kW]', power], \n",
    "                ['Battery Power [kW]', battery_power], \n",
    "                ['Energy Stored [MJ]', battery / 1000],\n",
    "                ['Total Regen [MJ]', total_regen / 1000],\n",
    "                ['Engine Power [kW]', engine_power], \n",
    "                ['Fuel Consumed [L]', fuel_consumed]])\n",
    "        \n",
    "    return {\n",
    "        'fuel_consumed': fuel_consumed, \n",
    "        'final_battery': battery, \n",
    "        'total_regen': total_regen / 1000,\n",
    "        'stop_times': stop_times,\n",
    "        'total_time': output[-1][0][1],\n",
    "        'total_dist': position,\n",
    "        'data': np.array(output)\n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def find_steady_state(data, margin=100):\n",
    "    initial_battery = 0\n",
    "    while True:\n",
    "        sim = run_simulation(data, hybrid=True, initial_battery=initial_battery)\n",
    "        battery = sim['final_battery']\n",
    "        if battery - margin <= initial_battery <= battery + margin:\n",
    "            return battery\n",
    "        initial_battery = battery"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def render_simulation(route, parameters, hybrid, initial_battery=0, time_limit=60, dt=0.05, reverse=False, winter=False):\n",
    "    sim = run_simulation(route, parameters, hybrid=hybrid, initial_battery=0, time_limit=60, dt=dt, reverse=reverse, winter=winter)\n",
    "    print('fuel rate: {:.2f} l/h'.format(sim['fuel_consumed'] / sim['total_time'] * 60))\n",
    "    kpl = sim['total_dist'] / 1000 / sim['fuel_consumed']\n",
    "    print('efficiency: {:.2f} km/l ({:.2f} mpg)'.format(kpl, kpl * 2.35))\n",
    "    print('distance: {:.2f} km'.format(sim['total_dist'] / 1000))\n",
    "    print('time: {:.2f}'.format(sim['total_time']))\n",
    "    print('ending battery: {:.0f} kJ'.format(sim['final_battery']))\n",
    "    print('total regen: {:.2f} MJ'.format(sim['total_regen']))\n",
    "\n",
    "\n",
    "    output = sim['data']\n",
    "    channels = output.shape[1] - 1\n",
    "    show_channels = [0, 1, 2, 3, 4, 5, 6, 7, 8]\n",
    "\n",
    "    tlim = (4.61, 6)\n",
    "\n",
    "    f, axes = pyplot.subplots(len(show_channels), figsize=(15, 3 * len(show_channels)))\n",
    "    for i, channel in enumerate(show_channels):\n",
    "        if len(show_channels) == 1:\n",
    "            axis = axes\n",
    "        else:\n",
    "            axis = axes[i]\n",
    "        axis.plot(output[int(tlim[0]*60/dt):int(tlim[1]*60/dt), 0, 1], \n",
    "                           output[int(tlim[0]*60/dt):int(tlim[1]*60/dt), channel + 1, 1])\n",
    "        axis.set(ylabel=output[0, channel + 1, 0])\n",
    "        for t in [t for t in sim['stop_times'] if tlim[0] < t < tlim[1]]:\n",
    "            axis.axvline(t, color='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "render_simulation(10, parameters, True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def parameter_partial(parameter_name, relative_epsilon=-.01):\n",
    "    alt_params = {name: parameter for name, parameter in parameters.items()}\n",
    "    alt_params[parameter_name] = parameters[parameter_name] * (1 + relative_epsilon)\n",
    "        \n",
    "    base_improvements = []\n",
    "    alt_improvements = []\n",
    "    \n",
    "    for route in ROUTES:\n",
    "        base_diesel = run_simulation(route, parameters, False, initial_battery=0, time_limit=60, dt=0.05, reverse=False, winter=False)\n",
    "        base_hybrid = run_simulation(route, parameters, True, initial_battery=0, time_limit=60, dt=0.05, reverse=False, winter=False)\n",
    "        alt_diesel = run_simulation(route, alt_params, False, initial_battery=0, time_limit=60, dt=0.05, reverse=False, winter=False)\n",
    "        alt_hybrid = run_simulation(route, alt_params, True, initial_battery=0, time_limit=60, dt=0.05, reverse=False, winter=False)\n",
    "        \n",
    "        base_diesel_rate = base_diesel['fuel_consumed'] / base_diesel['total_time']\n",
    "        base_hybrid_rate = base_hybrid['fuel_consumed'] / base_hybrid['total_time']\n",
    "        alt_diesel_rate = alt_diesel['fuel_consumed'] / alt_diesel['total_time']\n",
    "        alt_hybrid_rate = alt_hybrid['fuel_consumed'] / alt_hybrid['total_time']\n",
    "        \n",
    "        base_improvements.append(base_diesel_rate - base_hybrid_rate)\n",
    "        alt_improvements.append(alt_diesel_rate - alt_hybrid_rate)\n",
    "    \n",
    "#     print(base_improvements, alt_improvements)\n",
    "        \n",
    "    absolute_delta = 0\n",
    "    for bi, ai in zip(base_improvements, alt_improvements):\n",
    "        absolute_delta += abs((bi - ai) / bi)\n",
    "        \n",
    "#     print(absolute_delta)\n",
    "    \n",
    "    return absolute_delta / relative_epsilon / len(ROUTES)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for p in sorted(parameters.keys()):\n",
    "    print('{:<25} {:.3f}'.format(p, parameter_partial(p)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "winter = True\n",
    "\n",
    "for route in ROUTES:\n",
    "    diesel = run_simulation(route, parameters, False, 0, 60, 0.05, False, winter=winter)\n",
    "    h1 = run_simulation(route, parameters, True, 0, 60, 0.05, False, winter=winter)\n",
    "    hybrid = run_simulation(route, parameters, True, h1['final_battery'], 60, 0.05, False, winter=winter)\n",
    "    print('Route {}: Diesel: {:.2f} L/h, Hybrid: {:.2f} L/h'.format(\n",
    "            route, diesel['fuel_consumed'] / diesel['total_time'] * 60, hybrid['fuel_consumed'] / hybrid['total_time'] * 60))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "efficiencies = []\n",
    "stop_rates = []\n",
    "regen_rates = []\n",
    "\n",
    "for route in routes:\n",
    "    data_file = 'route_data/{}_data.json'.format(route)\n",
    "    data = json.load(open(data_file, 'r'))\n",
    "    diesel = run_simulation(data, hybrid=False, time_limit=120)\n",
    "    steady_battery = find_steady_state(data)\n",
    "    print(steady_battery)\n",
    "    hybrid = run_simulation(data, hybrid=True, initial_battery=steady_battery, time_limit=120)\n",
    "    stop_rate = len([p for p in data if p['stop']]) / hybrid['total_time'] * 60\n",
    "    diesel_rate = diesel['fuel_consumed'] / diesel['total_time'] * 60\n",
    "    hybrid_rate = hybrid['fuel_consumed'] / hybrid['total_time'] * 60\n",
    "    improvement = (diesel_rate - hybrid_rate) / diesel_rate * 100\n",
    "    regen_rate = hybrid['total_regen'] / hybrid['total_time'] / 60 * 1000\n",
    "    print('Route {}'.format(route))\n",
    "    print('diesel: {:.2f}, hybrid: {:.2f}, improvement: {:.2f}%, regen rate: {:.2f} kW, stops: {:.2f}/h'.format(\n",
    "            diesel_rate, hybrid_rate, improvement, regen_rate, stop_rate))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y = [6.38, 6.17, 6.12, 6.02, 5.67, 5.98]\n",
    "x = [5.14, 4.52, 4.14, 4.13, 3.98, 4.38]\n",
    "s = [89.0, 92.4, 100.0, 105.6, 89.4, 87.5]\n",
    "pyplot.bar(list(range(len(routes))), y, align='center')\n",
    "pyplot.xticks(list(range(len(routes))), [str(r) for r in routes])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# improvements = []\n",
    "# stop_rates = []\n",
    "avg_gradient = []\n",
    "\n",
    "for route in ROUTES:\n",
    "    print(route)\n",
    "    \n",
    "    data = load_route(route)\n",
    "    gradient = 0\n",
    "    for i in range(len(data) - 1):\n",
    "        gradient += abs(data[i]['gradient']) * (data[i + 1]['dist'] - data[i]['dist'])\n",
    "    avg_gradient.append(gradient / data[-1]['dist'])\n",
    "    \n",
    "#     stops = len([p for p in data if p['stop']])\n",
    "#     diesel = run_simulation(route, parameters, False, 0, 60, 0.05, False, False)\n",
    "#     h1 = run_simulation(route, parameters, True, 0, 60, 0.05, False, False)\n",
    "#     hybrid = run_simulation(route, parameters, True, h1['final_battery'], 60, 0.05, False, False)\n",
    "#     diesel_rate = diesel['fuel_consumed'] / diesel['total_time'] * 60\n",
    "#     hybrid_rate = hybrid['fuel_consumed'] / hybrid['total_time'] * 60\n",
    "#     improvements.append(diesel_rate - hybrid_rate)\n",
    "#     stop_rates.append(hybrid['total_time'] / stops)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pyplot.scatter(stop_rates, improvements)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.savetxt('out.csv', np.array(list(zip(improvements, stop_rates, avg_gradient))), delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "','.join(improvements)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "improvements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}