A workaround to divide a zip load into separate loads

PMD_zip_load_workaround.jl

using PowerModelsDistribution
using Ipopt
using Plots
const PMD = PowerModelsDistribution


NOTEBOOK_DIR = join(split(split(@__FILE__, "#==#")[1], "/")[1:end-1], "/")

data_eng = PMD.parse_file(
	"$NOTEBOOK_DIR/resources/lvtestcase_notrans.dss",
	transformations=[remove_all_bounds!]
);

PMD.reduce_lines!(data_eng);
data_eng["settings"]["sbase_default"] = 1.0*1E3/data_eng["settings"]["power_scale_factor"];
PMD.add_bus_absolute_vbounds!(
	data_eng,
	phase_lb_pu=0.9,
	phase_ub_pu=1.1,
	neutral_ub_pu=0.1
);

load_count = length(data_eng["load"])
load_ids = [id for (id, load) in data_eng["load"]]
for id in load_ids
    load = data_eng["load"][id]
    pd_nom = load["pd_nom"]
    qd_nom = load["qd_nom"]

    load["model"] = IMPEDANCE
    load["pd_cz"] = 0.2
    load["qd_cz"] = 0.2
    load["pd_nom_init"] = pd_nom
    load["qd_nom_init"] = qd_nom
    load["pd_nom"] = pd_nom * load["pd_cz"]
    load["qd_nom"] = qd_nom * load["qd_cz"]
    load["vm_nom"] = 0.22

    data_eng["load"]["$(id)_ci"] = deepcopy(load)
    data_eng["load"]["$(id)_ci"]["model"] = CURRENT
    data_eng["load"]["$(id)_ci"]["pd_ci"] = 0.3
    data_eng["load"]["$(id)_ci"]["qd_ci"] = 0.3
    data_eng["load"]["$(id)_ci"]["pd_nom"] = pd_nom * data_eng["load"]["$(id)_ci"]["pd_ci"]
    data_eng["load"]["$(id)_ci"]["qd_nom"] = qd_nom * data_eng["load"]["$(id)_ci"]["qd_ci"]

    data_eng["load"]["$(id)_cp"] = deepcopy(load)
    data_eng["load"]["$(id)_cp"]["model"] = POWER
    data_eng["load"]["$(id)_cp"]["pd_cp"] = 0.5
    data_eng["load"]["$(id)_cp"]["qd_cp"] = 0.5
    data_eng["load"]["$(id)_cp"]["pd_nom"] = pd_nom * data_eng["load"]["$(id)_cp"]["pd_cp"]
    data_eng["load"]["$(id)_cp"]["qd_nom"] = qd_nom * data_eng["load"]["$(id)_cp"]["qd_cp"]
end

res = solve_mc_opf(data_eng, ACPUPowerModel, Ipopt.Optimizer, setting=Dict("output"=>Dict("duals"=>true)));


##
pd_nom_init = [data_eng["load"]["$(id)"]["pd_nom_init"][1] for id in load_ids]
qd_nom_init = [data_eng["load"]["$(id)"]["qd_nom_init"][1] for id in load_ids]
pd_nom_cz = [data_eng["load"]["$(id)"]["pd_nom"][1] for id in load_ids]
qd_nom_cz = [data_eng["load"]["$(id)"]["qd_nom"][1] for id in load_ids]
pd_nom_ci = [data_eng["load"]["$(id)_ci"]["pd_nom"][1] for id in load_ids]
qd_nom_ci = [data_eng["load"]["$(id)_ci"]["qd_nom"][1] for id in load_ids]
pd_nom_cp = [data_eng["load"]["$(id)_cp"]["pd_nom"][1] for id in load_ids]
qd_nom_cp = [data_eng["load"]["$(id)_cp"]["qd_nom"][1] for id in load_ids]
pd_nom = pd_nom_cz .+ pd_nom_ci .+ pd_nom_cp
qd_nom = qd_nom_cz .+ qd_nom_ci .+ qd_nom_cp
vm_nom = 0.22

pd_cz = [res["solution"]["load"]["$(id)"]["pd"][1] for id in load_ids]
qd_cz = [res["solution"]["load"]["$(id)"]["qd"][1] for id in load_ids]
pd_ci = [res["solution"]["load"]["$(id)_ci"]["pd"][1] for id in load_ids]
qd_ci = [res["solution"]["load"]["$(id)_ci"]["qd"][1] for id in load_ids]
pd_cp = [res["solution"]["load"]["$(id)_cp"]["pd"][1] for id in load_ids]
qd_cp = [res["solution"]["load"]["$(id)_cp"]["qd"][1] for id in load_ids]
pd_tot = pd_cz .+ pd_ci .+ pd_cp
qd_tot = qd_cz .+ qd_ci .+ qd_cp

scatter(pd_tot, pd_nom, xlabel="solution", ylabel="input", label="pd")
scatter!(qd_tot, qd_nom, label="qd")

scatter(pd_cz./pd_nom_init, label="cz constant = 0.2", ylabel="constant factors")
scatter!(pd_ci./pd_nom_init, label="ci constant = 0.3")
scatter!(pd_cp./pd_nom_init, label="cp constant = 0.5")

connections = [data_eng["load"]["$id"]["connections"][1] for id in load_ids]
load_buses = [data_eng["load"]["$id"]["bus"] for id in load_ids]
vm_res = [res["solution"]["bus"][bus_id]["vm"][connections[i]] for (i,bus_id) in enumerate(load_buses)]

scatter(vm_res/vm_nom,  pd_cz ./ pd_nom/0.2, aspect_ratio=:equal, xlabel="vm p.u.", ylabel="pd p.u.", label="cz")
scatter!(vm_res/vm_nom, pd_ci ./ pd_nom/0.3, aspect_ratio=:equal, label="ci")
scatter!(vm_res/vm_nom, pd_cp ./ pd_nom/0.5, aspect_ratio=:equal, label="cp")