In [ ]:
from matplotlib import pyplot as plt
import pandas as pd
import sys
%matplotlib inline
sys.path.append("/home/docker/data/Optimization/")

In [ ]:
from district_opt import district_opt as dis

In [ ]:
df_case_3_A_energy = pd.read_csv("case_3_A_energy.csv", index_col = 0)
df_case_3_A_money = pd.read_csv("case_3_A_money.csv", index_col = 0)

df_case_3_B_energy = pd.read_csv("case_3_B_energy.csv", index_col = 0)
df_case_3_B_money = pd.read_csv("case_3_B_money.csv", index_col = 0)

df_case_3_C_energy = pd.read_csv("case_3_C_energy.csv", index_col = 0)
df_case_3_C_money = pd.read_csv("case_3_C_money.csv", index_col = 0)

In [ ]:
res_dis_M = dis.run_simulation(stop_time = 24*3600*6,
                             t_A = df_case_3_A_money.index, P_hvac_A = df_case_3_A_money["P_hvac"].values, P_batt_A = df_case_3_A_money["batt.P"].values,
                             t_B = df_case_3_B_money.index, P_hvac_B = df_case_3_B_money["P_hvac"].values, P_batt_B = df_case_3_B_money["batt.P"].values,
                             t_C = df_case_3_C_money.index, P_hvac_C = df_case_3_C_money["P_hvac"].values, P_batt_C = df_case_3_C_money["batt.P"].values)

In [ ]:
res_dis_E = dis.run_simulation(stop_time = 24*3600*6,
                             t_A = df_case_3_A_energy.index, P_hvac_A = df_case_3_A_energy["P_hvac"].values, P_batt_A = df_case_3_A_energy["batt.P"].values,
                             t_B = df_case_3_B_energy.index, P_hvac_B = df_case_3_B_energy["P_hvac"].values, P_batt_B = df_case_3_B_energy["batt.P"].values,
                             t_C = df_case_3_C_energy.index, P_hvac_C = df_case_3_C_energy["P_hvac"].values, P_batt_C = df_case_3_C_energy["batt.P"].values)

In [ ]:
def plot_temperatures(res):
    plt.plot(res["time"]/24/3600, res["buiA.Tmix"])
    plt.plot(res["time"]/24/3600, res["buiB.Tmix"])
    plt.plot(res["time"]/24/3600, res["buiC.Tmix"])
    plt.plot(res["time"][[0,-1]]/24/3600, [273.15+20,273.15+20], 'k')
    plt.plot(res["time"][[0,-1]]/24/3600, [273.15+24,273.15+24], 'k')
    plt.plot(res["time"]/24/3600, res["Tamb"])
    return

In [ ]:
plot_temperatures(res_dis_E)

In [ ]:
plot_temperatures(res_dis_M)

In [ ]:
def plot_solar_power(res):
    plt.plot(res["time"]/24/3600, res["buiA.pv.P"])
    plt.plot(res["time"]/24/3600, res["buiB.pv.P"])
    plt.plot(res["time"]/24/3600, res["buiC.pv.P"])

In [ ]:
plot_solar_power(res_dis_E)

In [ ]:
def plot_apparent_power(res):
    plt.plot(res["time"]/24/3600, res["buiA.S_pow"])
    plt.plot(res["time"]/24/3600, res["buiB.S_pow"])
    plt.plot(res["time"]/24/3600, res["buiC.S_pow"])
    return

In [ ]:
plot_apparent_power(res_dis_E)

In [ ]:
plot_apparent_power(res_dis_M)

In [ ]:
def plot_voltages(res):
    plt.plot(res["time"]/24/3600, res["buiA.Vrms"])
    plt.plot(res["time"]/24/3600, res["buiB.Vrms"])
    plt.plot(res["time"]/24/3600, res["buiC.Vrms"])
    plt.plot(res["time"][[0,-1]]/24/3600, [4800*0.95,4800*0.95], 'k')
    plt.plot(res["time"][[0,-1]]/24/3600, [4800*1,4800*1], 'k--')
    plt.plot(res["time"][[0,-1]]/24/3600, [4800*1.02,4800*1.02], 'k')
    return

In [ ]:
plot_voltages(res_dis_E)

In [ ]:
plot_voltages(res_dis_M)

In [ ]:
opt_res_dir = dis.run_optimization(res_dis_M, stop_time = 3600*24*6, opt_problem = 'ElectricNetwork.OptimizationDistrict_Money', n_e = 96)
#opt_res_dir = dis.run_optimization(opt_res_dir, n_e = 12)

In [ ]:
import matplotlib as mpl
mpl.rc('font',family='Times New Roman', size='13')

def plot_voltages_opt(res_E, res_M, res_opt_M, img_name = "aaa.png"):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(res_E["time"]/24/3600, res_E["buiA.Vrms"]/1e3,'b', alpha = 0.6, label="$V_{A,B,C}^E$")
    ax.plot(res_E["time"]/24/3600, res_E["buiB.Vrms"]/1e3,'b', alpha = 0.6)
    ax.plot(res_E["time"]/24/3600, res_E["buiC.Vrms"]/1e3,'b', alpha = 0.6)
    ax.plot(res_M["time"]/24/3600, res_M["buiA.Vrms"]/1e3,'r', alpha = 0.6, label="$V_{A,B,C}^M$")
    ax.plot(res_M["time"]/24/3600, res_M["buiB.Vrms"]/1e3,'r', alpha = 0.6)
    ax.plot(res_M["time"]/24/3600, res_M["buiC.Vrms"]/1e3,'r', alpha = 0.6)
    ax.plot(res_opt_M["time"]/24/3600, res_opt_M["buiA.Vrms"]/1e3,'g', alpha = 0.6, label="$V_{A,B,C}^{OPT}$")
    ax.plot(res_opt_M["time"]/24/3600, res_opt_M["buiB.Vrms"]/1e3,'g', alpha = 0.6)
    ax.plot(res_opt_M["time"]/24/3600, res_opt_M["buiC.Vrms"]/1e3,'g', alpha = 0.6)
    ax.plot(res_E["time"][[0,-1]]/24/3600, [4.800*1,4.800*1], 'k--',label="$V_n$")
    ax.fill_between(res_E["time"][[0,-1]]/24/3600, [4.800*1.02,4.800*1.02], [4.800*1.2,4.800*1.2], facecolor='red', alpha= 0.3, linewidth=0.0)
    plt.xlim([0, 5])
    plt.ylim([4.800*0.87, 4.800*1.12])
    ax.set_xlabel('Time [days]', size= '16')
    ax.set_ylabel('Voltage [kV]', size= '16')
    
    plt.tick_params(axis='both', which='major', labelsize=16)
    plt.tick_params(axis='both', which='minor', labelsize=16)
    plt.legend(loc="upper right", ncol=4)
    fig = plt.gcf()
    fig.set_size_inches(9,5)
    fig.savefig(img_name, dpi=300)
    plt.show()
    return

plot_voltages_opt(res_dis_E, res_dis_M, opt_res_dir, "voltage_comparison.png")

In [ ]:
print "Constraints on voltages","="*20
print opt_res_dir["buiA.Money"][-1]
print opt_res_dir["buiB.Money"][-1]
print opt_res_dir["buiC.Money"][-1]
print "No constraints on voltages",
print "="*20
print df_case_3_A_money["Money"].values[-1]
print df_case_3_B_money["Money"].values[-1]
print df_case_3_C_money["Money"].values[-1]
print "Difference","="*20
print (opt_res_dir["buiA.Money"][-1] - df_case_3_A_money["Money"].values[-1])/df_case_3_A_money["Money"].values[-1]
print (opt_res_dir["buiB.Money"][-1] - df_case_3_B_money["Money"].values[-1])/df_case_3_B_money["Money"].values[-1]
print (opt_res_dir["buiC.Money"][-1] - df_case_3_C_money["Money"].values[-1])/df_case_3_C_money["Money"].values[-1]

In [ ]:
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiA.Tmix"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiB.Tmix"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiC.Tmix"])
plt.plot(opt_res_dir["time"][[0,-1]]/24/3600, [273.15+20,273.15+20], 'k')
plt.plot(opt_res_dir["time"][[0,-1]]/24/3600, [273.15+24,273.15+24], 'k')
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["Tamb"])

In [ ]:
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiA.Vrms"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiB.Vrms"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiC.Vrms"])
plt.plot(opt_res_dir["time"][[0,-1]]/24/3600, [4800*0.95,4800*0.95], 'k')
plt.plot(opt_res_dir["time"][[0,-1]]/24/3600, [4800*1.05,4800*1.05], 'k')

In [ ]:
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiA.batt.SOC"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiB.batt.SOC"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiC.batt.SOC"])
plt.plot(opt_res_dir["time"][[0,-1]]/24/3600, [0.2,0.2], 'k')
plt.plot(opt_res_dir["time"][[0,-1]]/24/3600, [1.,1.], 'k')

In [ ]:
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiA.P_hvac"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiB.P_hvac"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiC.P_hvac"])

In [ ]:
opt_res_dir = dis.run_optimization(res_dis, "ElectricNetwork.OptimizationDistrict_Money", n_e = 12)

In [ ]:
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiA.batt.SOC"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiB.batt.SOC"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiC.batt.SOC"])
plt.plot(opt_res_dir["time"][[0,-1]]/24/3600, [0.2,0.2], 'k')
plt.plot(opt_res_dir["time"][[0,-1]]/24/3600, [1.,1.], 'k')

In [ ]:
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiA.P_hvac"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiB.P_hvac"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiC.P_hvac"])

In [ ]:
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiA.Vrms"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiB.Vrms"])
plt.plot(opt_res_dir["time"]/24/3600, opt_res_dir["buiC.Vrms"])
plt.plot(opt_res_dir["time"][[0,-1]]/24/3600, [4800*0.95,4800*0.95], 'k')
plt.plot(opt_res_dir["time"][[0,-1]]/24/3600, [4800*1.05,4800*1.05], 'k')

In [ ]:
plt.plot(opt_res_dir["buiA.P_el"])
plt.plot(opt_res_dir["buiA.pv.P"])
plt.plot(opt_res_dir["buiB.pv.P"])
plt.plot(opt_res_dir["buiC.pv.P"])

In [ ]: