# Test of the transmission and distribution combined model


## Import packages and choose a solver

using Printf
import PowerModels; const _PM = PowerModels
import PowerModelsACDC; const _PMACDC = PowerModelsACDC
import FlexPlan; const _FP = FlexPlan
import Cbc
optimizer = _FP.optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0)


# Include profile data from FlexPlan repository; you can of course also use your own data
include("../test/io/create_profile.jl")

## Input parameters

# Number of hourly optimization periods
number_of_hours = 12

# PowerModels, PowerModelsACDC and FlexPlan settings
s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => false)


## Transmission network instance (all data preparation except for `make_multinetwork` call)

t_file = "./test/data/case6/case6_2030.m" # Input case for transmission network

t_data = _FP.parse_file(t_file)
_FP.add_dimension!(t_data, :hour, number_of_hours)
_FP.add_dimension!(t_data, :scenario, Dict(1 => Dict{String,Any}("probability"=>1)), metadata = Dict{String,Any}("mc"=>true))
_FP.add_dimension!(t_data, :year, 1; metadata = Dict{String,Any}("scale_factor"=>1))
_FP.add_dimension!(t_data, :sub_nw, 1)
_FP.scale_data!(t_data)
t_data, t_loadprofile, t_genprofile = create_profile_data_italy!(t_data)
t_time_series = create_profile_data(number_of_hours, t_data, t_loadprofile, t_genprofile)


## Distribution network instance 1 (all data preparation except for `make_multinetwork` call)

d_file     = "test/data/cigre_mv_eu/cigre_mv_eu.m" # Input case for distribution networks
scale_load = 1.0 # Scaling factor of loads
scale_gen  = 1.0 # Scaling factor of generators

d_data_1 = _FP.parse_file(d_file)
_FP.add_dimension!(d_data_1, :hour, number_of_hours)
_FP.add_dimension!(d_data_1, :scenario, Dict(1 => Dict{String,Any}("probability"=>1)))
_FP.add_dimension!(d_data_1, :year, 1; metadata = Dict{String,Any}("scale_factor"=>1))
_FP.add_dimension!(d_data_1, :sub_nw, 1)
_FP.shift_ids!(d_data_1, _FP.dim_length(t_data))
_FP.scale_data!(d_data_1)
d_time_series = create_profile_data_cigre(d_data_1, number_of_hours; scale_load, scale_gen)
_FP.add_td_coupling_data!(t_data, d_data_1; t_bus = 1, sub_nw = 1) # The first distribution network is connected to bus 1 of transmission network.


## Distribution network instance 2 (all data preparation except for `make_multinetwork` call)

d_data_2 = deepcopy(d_data_1) # For simplicity, here a second distributon network is generated by duplicating the existing one
_FP.shift_ids!(d_data_2, _FP.dim_length(t_data))
_FP.add_td_coupling_data!(t_data, d_data_2; t_bus = 2, sub_nw = 1) # The second distribution network is connected to bus 2 of transmission network.


## Multinetwork data preparation

t_mn_data = _FP.make_multinetwork(t_data, t_time_series)
d_data_1 = _FP.make_multinetwork(d_data_1, d_time_series)
d_data_2 = _FP.make_multinetwork(d_data_2, d_time_series)
d_mn_data = _FP.merge_multinetworks!(d_data_1, d_data_2, :sub_nw)


## Solve problem

result = _FP.flex_tnep(t_mn_data, d_mn_data, _PM.DCPPowerModel, _FP.BFARadPowerModel, optimizer; setting=s)
@assert result["termination_status"] ∈ (_PM.OPTIMAL, _PM.LOCALLY_SOLVED) "$(result["optimizer"]) termination status: $(result["termination_status"])"


## Analyze result

t_nws = _FP.dim_length(t_mn_data)
d_nws = _FP.dim_length(d_mn_data)
sub_nws = _FP.dim_length(d_mn_data, :sub_nw)

printstyled("\n\n====================   Transmission network, first period   ====================\n\n", bold=true, color=:cyan)
_PM.print_summary(result["solution"]["nw"]["1"])

for sub_nw in 1:sub_nws
    printstyled("\n\n===================   Distribution network $sub_nw, first period   ===================\n\n", bold=true, color=:cyan)
    _PM.print_summary(result["solution"]["nw"]["$(sub_nw*t_nws+1)"])
end

printstyled("\n\n==========================   Power exchange at PCCs   ==========================\n\n", bold=true, color=:cyan)
println("Power in MW and MVar, positive if from transmission to distribution\n")
print("period ")
for sub_nw in 1:sub_nws
    @printf("%13s%10s", "p_dist$sub_nw", "q_dist$sub_nw")
end
println()
for t_nw in 1:t_nws
    @printf("%6i:", t_nw)
    for sub_nw in 1:sub_nws
        d_nw = sub_nw * t_nws + t_nw
        t_gen = _FP.dim_prop(d_mn_data, :sub_nw, sub_nw, "t_gen")
        t_res = result["solution"]["nw"]["$t_nw"]
        d_res = result["solution"]["nw"]["$d_nw"]
        t_p_in = t_res["gen"]["$t_gen"]["pg"] * t_res["baseMVA"]
        d_p_in = d_res["td_coupling"]["p"] * d_res["baseMVA"]
        @assert d_p_in ≈ -t_p_in
        d_q_in = d_res["td_coupling"]["q"] * d_res["baseMVA"]
        @printf("%13.3f%10.3f", d_p_in, d_q_in)
    end
    println()
end
