main.py

"""
Main OpenDSS Analysis Script
Modularized version of the original OpenDSS analysis
"""

import networkx as nx
import os
import opendssdirect as dss
import sys
import scipy.stats as st
import numpy as np
import math
import shutil
import pandas as pd
import matplotlib
from config import MATPLOTLIB_BACKEND
matplotlib.use(MATPLOTLIB_BACKEND)
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import matplotlib._color_data as mcd
import matplotlib.animation as ani
from scipy.stats import kurtosis, skew
from ditto.store import Store
from ditto.readers.json.read import Reader as JsonReader
from ditto.models import Line

# Import our modular components
from dss_checks import check_blown_fuses, check_voltage_violations, check_xfmr_overloads, check_line_overloads
from data_processing import process_timeseries_loads, process_peak_loads, get_all_feeders
from plotting_utils import do_plotting, create_voltage_histogram, create_voltage_percentiles
from output_utils import save_violation_results, save_summary_data, save_voltage_data
from timeseries_analysis import run_timeseries_analysis
from config import OUTPUT_FILES, ANALYSIS_FOLDER, LOAD_CURVES_FOLDER


def setup_folders(dataset, region, scenario, substation, feeder):
    """Setup output folders"""
    output_folder = os.path.join(dataset, region, 'scenarios', scenario, 'opendss', substation, feeder, ANALYSIS_FOLDER)
    if not os.path.isdir(output_folder):
        os.makedirs(output_folder)
    
    master_folder = os.path.join(dataset, region, 'scenarios', scenario, 'opendss', substation, feeder)
    
    return output_folder, master_folder


def run_opendss_analysis(dataset, region, year, scenario, substation, feeder, 
                        all_load_mappings, all_pv_mappings, profile_data, 
                        pv_profile_data, max_time, peak_loads):
    """Run OpenDSS analysis for a specific feeder"""
    
    print("substation:", substation, "feeder:", feeder)
    
    # Setup master file path
    master_file = os.path.join(dataset, region, 'scenarios', scenario, 'opendss_no_loadshapes', substation, feeder, 'Master.dss')
    if year == 'peak':
        master_file = os.path.join(dataset, region, 'scenarios', scenario, 'opendss', substation, feeder, 'Master.dss')

    print('starting initial load for', dataset, region, substation, feeder, scenario)
    result = dss.run_command("Redirect " + master_file)
    print('finished initial load')
    print(result, flush=True)
    
    # Setup folders
    output_folder, master_folder = setup_folders(dataset, region, scenario, substation, feeder)
    
    # Initialize load values
    max_total = 0
    max_res = 0
    max_com = 0
    all_kw = []
    all_res = []
    all_com = []

    if year == 'peak':
        max_total, max_res, max_com = peak_loads[(substation, feeder)]
    
    # Process loads for timeseries analysis
    if not year == 'peak':
        dss.Loads.First()
        while True:
            name = dss.Loads.Name()
            name_base = name
            multiplier = 1
            if name.endswith('_1') or name.endswith('_2'):  # center tap loads
                name_base = name[:-2]
                multiplier = 0.5
        
            parquet_name = all_load_mappings[name][0]
            input_mult = all_load_mappings[name][1]
            if input_mult != multiplier:
                print(f'{parquet_name} has mult inconsistencies')
            if 'mesh' in parquet_name:
                parquet_name = 'mesh'
            kw = profile_data[parquet_name][0].iloc[max_time] * multiplier
            kvar = profile_data[parquet_name][1].iloc[max_time] * multiplier
        
            max_total += kw
            all_kw.append(profile_data[parquet_name][0] * multiplier)
            if parquet_name.startswith('res_'):
                max_res += kw
                all_res.append(profile_data[parquet_name][0] * multiplier)

            if parquet_name.startswith('com_') or parquet_name == 'mesh':
                max_com += kw
                all_com.append(profile_data[parquet_name][0] * multiplier)

            res1 = dss.Loads.kW(kw)
            res2 = dss.Loads.kvar(kvar)
        
            if not dss.Loads.Next() > 0:
                break
    
    # Process PV systems
    if not year == 'peak':
        first_pv = dss.PVsystems.First()
        while True and first_pv > 0:
            name = dss.PVsystems.Name()
            csv_name = all_pv_mappings[name]
            irradiance = pv_profile_data[csv_name].iloc[max_time]
            res1 = dss.PVsystems.Irradiance(irradiance/1000)
        
            if not dss.PVsystems.Next() > 0:
                break
    
    # Solve the system
    if not year == 'peak':
        print('starting opendss run')
        result = dss.run_command(f"Solve")
        print('finished opendss run')
        print(result)
    
    # Calculate losses and power
    all_losses = dss.Circuit.Losses()[0]/1000
    total_power = dss.Circuit.TotalPower()[0]*-1
    if total_power == 0:
        percent_losses = 0
    else:
        percent_losses = 100*all_losses/total_power
    if total_power < 0:
        percent_losses = percent_losses*-1
    
    # Check system violations
    undervoltages_dict, overvoltages_dict = check_voltage_violations()
    xfmr_overloads_dict, unloaded_xfmrs_dict = check_xfmr_overloads()
    line_overloads_dict, unloaded_lines_dict = check_line_overloads()
    blown_fuses = check_blown_fuses()
    
    # Save violation results
    save_violation_results(output_folder, undervoltages_dict, overvoltages_dict, 
                          xfmr_overloads_dict, line_overloads_dict, blown_fuses)
    
    # Save summary data (need peak time values from processing)
    peak_total_day = ''
    peak_total_hour = '' 
    peak_total_min = ''
    
    save_summary_data(output_folder, all_losses, total_power, percent_losses, 
                     max_total, max_com, max_res, peak_total_day, peak_total_hour, peak_total_min)
    
    print(f'All losses: {all_losses}')
    print(f'Total power {total_power}')
    print(f'Percentage Losses {percent_losses}')
    
    # Calculate and save voltage data
    global_max = -1000
    global_min = 100000
    names = dss.Circuit.AllBusNames()
    dss_vals_avg = []
    sorted_dss_vals_avg = []
    deenergized_switch_nodes = []
    
    for name in names:
        dss.Circuit.SetActiveBus(name)
        dss_pus = dss.Bus.PuVoltage()
        tot = 0
        cnt = 0
        for i in range(int(len(dss_pus)/2)):
            mag = abs(complex(dss_pus[2*i], dss_pus[2*i+1]))
            if mag > 0:
                tot += mag
                cnt += 1
        if cnt == 0:
            tot = 0
        else:
            tot = tot/cnt
        if tot > 0:
            sorted_dss_vals_avg.append((tot, name))
            dss_vals_avg.append(tot)
    
            if tot < global_min:
                global_min = tot
            if tot > global_max:
                global_max = tot
        else:
            deenergized_switch_nodes.append(name)
    
    # Save voltage data
    save_voltage_data(output_folder, sorted_dss_vals_avg, deenergized_switch_nodes, master_folder)
    
    try:
        sorted_dss_vals_avg_df = pd.DataFrame(sorted(sorted_dss_vals_avg))
        if len(sorted_dss_vals_avg_df) > 0:
            print(f'Min Voltage: {sorted_dss_vals_avg_df[0].iloc[0]}')
            print(f'Max Voltage: {sorted_dss_vals_avg_df[0].iloc[-1]}')
    except:
        print("Could not display voltage range")
    
    # Create plots
    if len(dss_vals_avg) > 0:
        create_voltage_histogram(dss_vals_avg, global_min, global_max, region, substation, feeder, output_folder)
        create_voltage_percentiles(dss_vals_avg, region, substation, feeder, output_folder)
    
    return all_kw, all_res, all_com


def main():
    """Main function"""
    if len(sys.argv) != 7:
        print("Usage: python main.py <dataset> <region> <year> <scenario> <run_all> <delete_folders>")
        sys.exit(1)
    
    dataset = sys.argv[1]
    region = sys.argv[2]
    year = sys.argv[3]
    scenario = sys.argv[4]
    run_all = sys.argv[5]
    delete_folders = sys.argv[6]

    raw_dataset = dataset
    dataset = os.path.join(year, dataset)
    
    # Check if master file exists
    master_file = os.path.join(dataset, region, 'scenarios', scenario, 'opendss_no_loadshapes', 'Master.dss')
    if year == 'peak':
        master_file = os.path.join(dataset, region, 'scenarios', scenario, 'opendss', 'Master.dss')
    if not os.path.exists(master_file):
        print(f'{region}, {dataset} is missing. skipping...')
        return

    # Process data based on year type
    if year == 'peak':
        peak_loads = process_peak_loads(dataset, region, scenario)
        all_load_mappings = {}
        all_pv_mappings = {}
        profile_data = {}
        pv_profile_data = {}
        max_time = 0
        total_load = None
    else:
        (all_load_mappings, all_pv_mappings, profile_data, pv_profile_data, 
         total_load, max_time, peak_total_day, peak_total_hour, peak_total_min,
         region_all_kw, all_kvar) = process_timeseries_loads(dataset, region, scenario)
        peak_loads = {}

    # Get all feeders to process
    base_folder = os.path.join(dataset, region, 'scenarios', scenario, 'opendss_no_loadshapes')
    all_feeders = get_all_feeders(base_folder)

    # Process each feeder
    for substation, feeder in all_feeders:
        try:
            all_kw, all_res, all_com = run_opendss_analysis(
                dataset, region, year, scenario, substation, feeder,
                all_load_mappings, all_pv_mappings, profile_data, pv_profile_data,
                max_time, peak_loads
            )
            
            # Run timeseries analysis if requested
            if run_all == 'timeseries' and year != 'peak':
                output_folder = os.path.join(dataset, region, 'scenarios', scenario, 'opendss', 
                                           substation, feeder, ANALYSIS_FOLDER)
                run_timeseries_analysis(output_folder, all_load_mappings, all_pv_mappings, 
                                       profile_data, pv_profile_data, year)
            
        except Exception as e:
            print(f"Error processing {substation}/{feeder}: {e}")
            continue

    # Clean up folders if requested
    if delete_folders == 'Delete':
        delete_base = os.path.join(dataset, region, 'scenarios', scenario)
        delete_json = os.path.join(delete_base, 'json_opendss')
        delete_prelim = os.path.join(delete_base, 'opendss_prelim')
        if os.path.exists(delete_json):
            print(f'Deleting folder {delete_json}')
            shutil.rmtree(delete_json)
        if os.path.exists(delete_prelim):
            print(f'Deleting folder {delete_prelim}')
            shutil.rmtree(delete_prelim)
    else:
        print('Not deleting any folders')
    
    # Generate load curves if this is base_timeseries scenario
    if scenario == 'base_timeseries' and year != 'peak' and total_load is not None:
        for substation, feeder in all_feeders:
            load_curve_subfolder = region
            if substation != '':
                load_curve_subfolder = load_curve_subfolder + '__' + substation
            if feeder != '':
                load_curve_subfolder = load_curve_subfolder + '__' + feeder
            load_curve_folder = os.path.join(dataset, region, LOAD_CURVES_FOLDER, load_curve_subfolder)
            if not os.path.isdir(load_curve_folder):
                os.makedirs(load_curve_folder)
        
            num_days = 365
            print('Writing daily plots')
            # Sequential processing for plotting
            total_res = sum(all_res) if all_res else pd.Series([0] * len(total_load))
            total_com = sum(all_com) if all_com else pd.Series([0] * len(total_load))
            
            for day in range(num_days):
                do_plotting(total_load, total_res, total_com, region, substation, feeder, day, load_curve_folder)


if __name__ == '__main__':
    main()