From 35f349439a8bbc177b1e639715c6ff7b5c767f15 Mon Sep 17 00:00:00 2001
From: jtabarez <53569964+jtabarez@users.noreply.github.com>
Date: Mon, 11 Mar 2024 14:30:52 -0600
Subject: [PATCH] gfli test added

---
 src/core/admittance_matrix.jl |  19 ++-
 src/data_model/components.jl  |   2 +-
 src/data_model/eng2math.jl    |   2 +-
 src/io/dss/dss2eng.jl         |  19 +++
 src/io/opendss.jl             |   2 +-
 test/data/dist/gfmi_test.dss  |  58 +++++++
 test/gfli.jl                  | 294 ++++++++++++++++++++++++++++++++++
 7 files changed, 391 insertions(+), 5 deletions(-)
 create mode 100644 test/data/dist/gfmi_test.dss
 create mode 100644 test/gfli.jl

diff --git a/src/core/admittance_matrix.jl b/src/core/admittance_matrix.jl
index f715522..929eeff 100644
--- a/src/core/admittance_matrix.jl
+++ b/src/core/admittance_matrix.jl
@@ -433,7 +433,7 @@ end
 function build_mc_delta_current_generator!(delta_i, v, data)
     for (_, gen) in data["gen"]
         if occursin("generator", gen["source_id"])
-            if gen["gen_model"] == 1
+            if gen["gen_model"] == 1 
                 calc_delta_current_generator!(gen, delta_i, v, data)
             end
         end
@@ -494,7 +494,6 @@ function build_mc_delta_current_inverter!(delta_i, v, data, z_matrix)
         if occursin("solar.", gen["source_id"])
             if gen["pv_model"] == 1
                 if gen["grid_forming"]
-                    println(oooooo)
                     calc_mc_delta_current_gfmi!(gen, delta_i, v, data)
                 else
                     calc_mc_delta_current_gfli!(gen, delta_i, v, data)
@@ -595,6 +594,7 @@ function update_mc_delta_current_vector(model, v)
     delta_i = zeros(Complex{Float64}, n, 1)
     update_mc_delta_current_generator!(delta_i, v, model.data)
     update_mc_delta_current_load!(delta_i, v, model.data)
+    update_mc_delta_current_inverter!(delta_i, v, model.data)
     return _SP.sparse(delta_i)
 end
 
@@ -618,4 +618,19 @@ function update_mc_delta_current_load!(delta_i, v, data)
             end
         end
     end
+end
+
+
+function update_mc_delta_current_inverter!(delta_i, v, data)
+    for (_, gen) in data["gen"]
+        if occursin("solar.", gen["source_id"])
+            if gen["pv_model"] == 1
+                if gen["grid_forming"]
+                    calc_mc_delta_current_gfmi!(gen, delta_i, v, data)
+                else
+                    calc_mc_delta_current_gfli!(gen, delta_i, v, data)
+                end
+            end
+        end
+    end
 end
\ No newline at end of file
diff --git a/src/data_model/components.jl b/src/data_model/components.jl
index a44dfd2..76b7836 100644
--- a/src/data_model/components.jl
+++ b/src/data_model/components.jl
@@ -495,7 +495,7 @@ function _map_eng2math_mc_admittance_2w_transformer!(transformer::Dict{String,<:
                     end
                 end
                 if w == 1
-                    if transformer["leadlag"] == "lag"
+                    if transformer["leadlag"] == "lead"
                         a[1,1] = a[1,6] = a[2,5] = a[2,10] = a[3,9] = a[3,2] = 1
                     else
                         a[1,1] = a[1,10] = a[2,2] = a[2,5] = a[3,6] = a[3,9] = 1
diff --git a/src/data_model/eng2math.jl b/src/data_model/eng2math.jl
index 1278798..dcb5324 100644
--- a/src/data_model/eng2math.jl
+++ b/src/data_model/eng2math.jl
@@ -157,7 +157,7 @@ end
 
 "field/values to passthrough from the ENGINEERING to MATHEMATICAL data models"
 const _pmp_eng2math_passthrough = Dict{String,Vector{String}}(
-        "generator" => String["zr", "zx", "grid_forming", "gen_model", "xdp", "rp", "xdpp", "vnom_kv"],
+        "generator" => String["zr", "zx", "gen_model", "xdp", "rp", "xdpp", "vnom_kv"],
         "solar" => String["i_max", "solar_max", "kva", "pf", "grid_forming", "balanced", "vminpu", "transformer", "type", "pv_model", "transformer_id"],
         "voltage_source" => String["zr", "zx"],
         "load" => String["vminpu", "vmaxpu"],
diff --git a/src/io/dss/dss2eng.jl b/src/io/dss/dss2eng.jl
index 708edfb..4961a3f 100644
--- a/src/io/dss/dss2eng.jl
+++ b/src/io/dss/dss2eng.jl
@@ -459,3 +459,22 @@ function _dss2eng_solar_transformer(data_eng::Dict{String,<:Any}, data_dss::Dict
         end
     end
 end
+
+
+function _dss2eng_issues!(data_eng::Dict{String,<:Any}, data_dss::Dict{String,<:Any})
+#     println(keys(data_eng))
+#     # println(data_eng["bus"])
+#     if haskey(data_eng, "line")
+#         for (id, line) in data_eng["line"]
+#             if haskey(line["dss"], "geometry")
+#                 # println(line)
+#                 # calc_line_paramters!(line, data_dss)
+#                 if length(line["rs"]) != length(line["f_connections"])^2
+
+#                     println(line)
+#                 end
+#             end
+#         end
+#     end
+#     println(oooo)
+end
diff --git a/src/io/opendss.jl b/src/io/opendss.jl
index 4c1d776..efe5a79 100644
--- a/src/io/opendss.jl
+++ b/src/io/opendss.jl
@@ -9,7 +9,7 @@ data model adding some extra data transformations specific to fault studies.
 `kwargs` can be any valid keyword arguments from PowerModelsDistribution's `parse_file` function.
 """
 function parse_opendss(file::String; method::Union{String,Missing}=missing, kwargs...)
-    pm_data = _PMD.parse_file(file; dss2eng_extensions=[_dss2eng_solar_dynamics!, _dss2eng_load_dynamics!, _dss2eng_gen_dynamics!, _dss2eng_transformer_dynamics!, _dss2eng_transformers!], kwargs...)
+    pm_data = _PMD.parse_file(file; dss2eng_extensions=[_dss2eng_solar_dynamics!, _dss2eng_load_dynamics!, _dss2eng_gen_dynamics!, _dss2eng_transformer_dynamics!, _dss2eng_transformers!, _dss2eng_issues!], kwargs...)
     if ismissing(method)
         pm_data["method"] = "PMD"
     else
diff --git a/test/data/dist/gfmi_test.dss b/test/data/dist/gfmi_test.dss
new file mode 100644
index 0000000..4abe9c9
--- /dev/null
+++ b/test/data/dist/gfmi_test.dss
@@ -0,0 +1,58 @@
+Clear
+New Circuit.3Bus_example
+!  define a really stiff source
+~ basekv=4.16 BaseMVA=.5  pu=0.9959  r1=.7  x1=3.5
+
+!Define Linecodes
+
+
+New linecode.556MCM nphases=3 basefreq=60  ! ohms per 5 mile
+~ rmatrix = ( 0.1000 | 0.0400    0.1000 |  0.0400    0.0400    0.1000)
+~ xmatrix = ( 0.0583 |  0.0233    0.0583 | 0.0233    0.0233    0.0583)
+~ cmatrix = (50.92958178940651  | -0  50.92958178940651 | -0 -0 50.92958178940651  ) ! small capacitance
+
+
+New linecode.4/0QUAD nphases=3 basefreq=60  ! ohms per 100ft
+~ rmatrix = ( 0.1167 | 0.0467    0.1167 | 0.0467    0.0467    0.1167)
+~ xmatrix = (0.0667  |  0.0267    0.0667  |  0.0267    0.0267    0.0667 )
+~ cmatrix = (50.92958178940651  | -0  50.92958178940651 | -0 -0 50.92958178940651  )  ! small capacitance
+
+!Define lines
+
+New Line.OHLine  bus1=sourcebus.1.2.3  Primary.1.2.3  linecode = 556MCM   length=1  ! 5 mile line
+New Line.Quad    Bus1=Primary.1.2.3  loadbus.1.2.3  linecode = 4/0QUAD  length=1   ! 100 ft
+New Line.PV_line    Bus1=PV_bus.1.2.3 Primary.1.2.3  linecode = 4/0QUAD  length=1   ! 100 ft
+
+!Loads - single phase
+
+!New Load.L1 phases=3  bus1=loadbus.1.2.3  Conn=Wye kW=160  kvar=130  model=1 kV=4.16
+New Load.L1 phases=3  bus1=loadbus.1.2.3  Conn=Wye kW=18  kvar=8  model=1 kV=4.16
+
+New Transformer.TX1 windings=2 phases=3 Buses=[PV_bus1 PV_bus]
+~ Conns=[Delta Wye]
+~ kVs=[.4 4.16]
+~ kVAs=[25 25]
+~ %Rs=[1 2]
+~ xhl=5
+~ %noloadloss=0.0
+~ %imag=0.0
+~ leadlag=lag
+
+New PVSystem.PV1 phases=3 bus1=PV_bus1.1.2.3 kv=0.4 conn=wye irradiance=1 Pmpp=25 kva=25 pf=.95 balanced=true limitcurrent=true vminpu=0.9090909
+
+
+Set voltagebases=[4.16 0.4]
+Set tolerance=0.000001
+set defaultbasefreq=60
+Calcvoltagebases
+solve
+
+!batchedit load..* enable=false
+!batchedit capacitor..* enabled=false
+new fault.test phases=3 bus1=loadbus.1.2.3 r=.0001 enabled=true
+
+Solve mode=snap
+
+set mode=dynamic controlmode=time
+set stepsize=.0002s number=800
+solve
diff --git a/test/gfli.jl b/test/gfli.jl
new file mode 100644
index 0000000..4a8c8d6
--- /dev/null
+++ b/test/gfli.jl
@@ -0,0 +1,294 @@
+
+
+using PowerModelsProtection
+
+import PowerModels
+import PowerModelsDistribution
+
+import Ipopt
+
+using Printf
+
+const PMD = PowerModelsDistribution
+
+PowerModelsDistribution.silence!()
+
+ipopt_solver = optimizer_with_attributes(Ipopt.Optimizer, "tol"=>1e-6, "print_level"=>0)
+
+data = parse_file("../test/data/dist/gfmi_test.dss")
+# Dict{Any, Int64} with 5 entries:
+#   "primary"   => 1
+#   "sourcebus" => 2
+#   "loadbus"   => 3
+#   "pv_bus"    => 4
+#   "pv_bus1"   => 5
+
+data["solar"]["pv1"]["transformer_id"] = "tx1"
+data["solar"]["pv1"]["grid_forming"] = false
+data["solar"]["pv1"]["pv_model"] = 1
+
+data_math = transform_admittance_data_model(data)
+
+model = instantiate_mc_admittance_model(data_math)        
+
+"""
+    Normal power flow
+"""
+y = deepcopy(model.y)
+
+v = compute_mc_pf(y, model)
+data = model.data
+gen = data["gen"]["1"]
+transformer = data["transformer"][gen["transformer_id"]]
+f_bus = data["bus"]["$(transformer["f_bus"])"]
+t_bus = data["bus"]["$(transformer["t_bus"])"]
+y = transformer["p_matrix"][5:8,1:8]
+_v = zeros(Complex{Float64}, 8, 1)
+indx = 1
+for (_j, j) in enumerate(f_bus["terminals"])
+    if haskey(data["admittance_map"], (f_bus["bus_i"], j))
+        _v[indx, 1] = v[data["admittance_map"][(f_bus["bus_i"], j)], 1]
+    else
+        _v[indx, 1] = 0.0
+    end
+    global indx += 1
+end
+for (_j, j) in enumerate(t_bus["terminals"])
+    if haskey(data["admittance_map"], (t_bus["bus_i"], j))
+        _v[indx, 1] = v[data["admittance_map"][(t_bus["bus_i"], j)], 1]
+    else
+        _v[indx, 1] = 0.0
+    end
+    global indx += 1
+end
+i_abc = (transformer["p_matrix"][1:4,1:8]*_v)
+i_line = (transformer["p_matrix"][5:8,1:8]*_v)
+i_012 = inv(PowerModelsProtection._A)*i_abc[1:3]
+v_012 = inv(PowerModelsProtection._A)*_v[1:3]
+println("_____power flow____")
+@printf("inverter phase currents (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(i_abc[1]), angle(i_abc[1])*180/pi, abs(i_abc[2]), angle(i_abc[2])*180/pi,abs(i_abc[3]), angle(i_abc[3])*180/pi)
+@printf("inverter phase voltages (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(_v[1]), angle(_v[1])*180/pi, abs(_v[2]), angle(_v[2])*180/pi, abs(_v[3]), angle(_v[3])*180/pi, abs(_v[4]), angle(_v[4])*180/pi)
+@printf("inverter sequence currents (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(i_012[1]), angle(i_012[1])*180/pi, abs(i_012[2]), angle(i_012[2])*180/pi,abs(i_012[3]), angle(i_012[3])*180/pi)
+@printf("inverter sequence voltages (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(v_012[1]), angle(v_012[1])*180/pi, abs(v_012[2]), angle(v_012[2])*180/pi,abs(v_012[3]), angle(v_012[3])*180/pi)
+
+i_012 = inv(PowerModelsProtection._A)*i_line[1:3]
+v_012 = inv(PowerModelsProtection._A)*_v[5:7]
+@printf("secondary side transformer phase currents (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(i_line[1]), angle(i_line[1])*180/pi, abs(i_line[2]), angle(i_line[2])*180/pi, abs(i_line[3]), angle(i_line[3])*180/pi, abs(i_line[4]), angle(i_line[4])*180/pi)
+@printf("secondary side transformer phase voltages (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(_v[5]), angle(_v[5])*180/pi, abs(_v[6]), angle(_v[6])*180/pi, abs(_v[7]), angle(_v[7])*180/pi, abs(_v[8]), angle(_v[8])*180/pi)
+@printf("secondary side transformer sequence currents (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(i_012[1]), angle(i_012[1])*180/pi, abs(i_012[2]), angle(i_012[2])*180/pi,abs(i_012[3]), angle(i_012[3])*180/pi)
+@printf("secondary side transformer sequence voltages (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(v_012[1]), angle(v_012[1])*180/pi, abs(v_012[2]), angle(v_012[2])*180/pi,abs(v_012[3]), angle(v_012[3])*180/pi)
+
+phase_resistance = .0001
+ground_resistance = .0001
+
+
+"""
+    ll- fault at load bus
+"""
+
+gp = 1 / phase_resistance
+gf = 1 / ground_resistance
+Gf_ll = zeros(Real, 2, 2)
+
+for i in 1:2
+    for j in 1:2
+        if i == j
+            Gf_ll[i,j] = 1 / phase_resistance
+        else
+            Gf_ll[i,j] = -1 / phase_resistance
+        end
+    end
+end
+
+Gf_ll[1,1] += gf
+Gf_ll[2,2] += gf
+
+                                                           
+bus = model.data["bus"]["3"]
+y = deepcopy(model.y)
+for i_indx in 1:2
+    for j_indx in 1:2
+        i = model.data["admittance_map"][(bus["bus_i"],bus["terminals"][i_indx])]
+        j = model.data["admittance_map"][(bus["bus_i"],bus["terminals"][j_indx])]
+        y[i,j] += Gf_ll[i_indx,j_indx]
+    end
+end
+
+v = compute_mc_pf(y, model)
+data = model.data
+gen = data["gen"]["1"]
+transformer = data["transformer"][gen["transformer_id"]]
+f_bus = data["bus"]["$(transformer["f_bus"])"]
+t_bus = data["bus"]["$(transformer["t_bus"])"]
+y = transformer["p_matrix"][5:8,1:8]
+_v = zeros(Complex{Float64}, 8, 1)
+indx = 1
+for (_j, j) in enumerate(f_bus["terminals"])
+    if haskey(data["admittance_map"], (f_bus["bus_i"], j))
+        _v[indx, 1] = v[data["admittance_map"][(f_bus["bus_i"], j)], 1]
+    else
+        _v[indx, 1] = 0.0
+    end
+    global indx += 1
+end
+for (_j, j) in enumerate(t_bus["terminals"])
+    if haskey(data["admittance_map"], (t_bus["bus_i"], j))
+        _v[indx, 1] = v[data["admittance_map"][(t_bus["bus_i"], j)], 1]
+    else
+        _v[indx, 1] = 0.0
+    end
+    global indx += 1
+end
+i_abc = (transformer["p_matrix"][1:4,1:8]*_v)
+i_line = (transformer["p_matrix"][5:8,1:8]*_v)
+i_012 = inv(PowerModelsProtection._A)*i_abc[1:3]
+v_012 = inv(PowerModelsProtection._A)*_v[1:3]
+println("_____ll-fault____")
+@printf("inverter phase currents (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(i_abc[1]), angle(i_abc[1])*180/pi, abs(i_abc[2]), angle(i_abc[2])*180/pi,abs(i_abc[3]), angle(i_abc[3])*180/pi)
+@printf("inverter phase voltages (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(_v[1]), angle(_v[1])*180/pi, abs(_v[2]), angle(_v[2])*180/pi, abs(_v[3]), angle(_v[3])*180/pi, abs(_v[4]), angle(_v[4])*180/pi)
+@printf("inverter sequence currents (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(i_012[1]), angle(i_012[1])*180/pi, abs(i_012[2]), angle(i_012[2])*180/pi,abs(i_012[3]), angle(i_012[3])*180/pi)
+@printf("inverter sequence voltages (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(v_012[1]), angle(v_012[1])*180/pi, abs(v_012[2]), angle(v_012[2])*180/pi,abs(v_012[3]), angle(v_012[3])*180/pi)
+
+i_012 = inv(PowerModelsProtection._A)*i_line[1:3]
+v_012 = inv(PowerModelsProtection._A)*_v[5:7]
+@printf("secondary side transformer phase currents (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(i_line[1]), angle(i_line[1])*180/pi, abs(i_line[2]), angle(i_line[2])*180/pi, abs(i_line[3]), angle(i_line[3])*180/pi, abs(i_line[4]), angle(i_line[4])*180/pi)
+@printf("secondary side transformer phase voltages (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(_v[5]), angle(_v[5])*180/pi, abs(_v[6]), angle(_v[6])*180/pi, abs(_v[7]), angle(_v[7])*180/pi, abs(_v[8]), angle(_v[8])*180/pi)
+@printf("secondary side transformer sequence currents (0, 1 ,2) = %f @ %F, %f @ %f, %f @ %F\n", abs(i_012[1]), angle(i_012[1])*180/pi, abs(i_012[2]), angle(i_012[2])*180/pi,abs(i_012[3]), angle(i_012[3])*180/pi)
+@printf("secondary side transformer sequence voltages (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(v_012[1]), angle(v_012[1])*180/pi, abs(v_012[2]), angle(v_012[2])*180/pi,abs(v_012[3]), angle(v_012[3])*180/pi)
+
+
+
+"""
+    3pg- fault at load bus
+"""
+
+Gf_3p = zeros(Real, 3, 3)
+for i in 1:3
+    for j in 1:3
+        if i != j
+            Gf_3p[i,j] = -1/phase_resistance
+        else
+            Gf_3p[i,j] = 2 * (1/phase_resistance)
+        end
+    end
+end
+                                                          
+bus = model.data["bus"]["3"]
+indx = (1,2)
+y = deepcopy(model.y)
+for i_indx in 1:3
+    for j_indx in 1:3
+        i = model.data["admittance_map"][(bus["bus_i"],bus["terminals"][i_indx])]
+        j = model.data["admittance_map"][(bus["bus_i"],bus["terminals"][j_indx])]
+        y[i,j] += Gf_3p[i_indx,j_indx]
+    end
+end
+
+v = compute_mc_pf(y, model)
+data = model.data
+gen = data["gen"]["1"]
+transformer = data["transformer"][gen["transformer_id"]]
+f_bus = data["bus"]["$(transformer["f_bus"])"]
+t_bus = data["bus"]["$(transformer["t_bus"])"]
+y = transformer["p_matrix"][5:8,1:8]
+_v = zeros(Complex{Float64}, 8, 1)
+indx = 1
+for (_j, j) in enumerate(f_bus["terminals"])
+    if haskey(data["admittance_map"], (f_bus["bus_i"], j))
+        _v[indx, 1] = v[data["admittance_map"][(f_bus["bus_i"], j)], 1]
+    else
+        _v[indx, 1] = 0.0
+    end
+    global indx += 1
+end
+for (_j, j) in enumerate(t_bus["terminals"])
+    if haskey(data["admittance_map"], (t_bus["bus_i"], j))
+        _v[indx, 1] = v[data["admittance_map"][(t_bus["bus_i"], j)], 1]
+    else
+        _v[indx, 1] = 0.0
+    end
+    global indx += 1
+end
+i_abc = (transformer["p_matrix"][1:4,1:8]*_v)
+i_line = (transformer["p_matrix"][5:8,1:8]*_v)
+i_012 = inv(PowerModelsProtection._A)*i_abc[1:3]
+v_012 = inv(PowerModelsProtection._A)*_v[1:3]
+println("_____3pg-fault____")
+@printf("inverter phase currents (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(i_abc[1]), angle(i_abc[1])*180/pi, abs(i_abc[2]), angle(i_abc[2])*180/pi,abs(i_abc[3]), angle(i_abc[3])*180/pi)
+@printf("inverter phase voltages (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(_v[1]), angle(_v[1])*180/pi, abs(_v[2]), angle(_v[2])*180/pi, abs(_v[3]), angle(_v[3])*180/pi, abs(_v[4]), angle(_v[4])*180/pi)
+@printf("inverter sequence currents (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(i_012[1]), angle(i_012[1])*180/pi, abs(i_012[2]), angle(i_012[2])*180/pi,abs(i_012[3]), angle(i_012[3])*180/pi)
+@printf("inverter sequence voltages (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(v_012[1]), angle(v_012[1])*180/pi, abs(v_012[2]), angle(v_012[2])*180/pi,abs(v_012[3]), angle(v_012[3])*180/pi)
+
+i_012 = inv(PowerModelsProtection._A)*i_line[1:3]
+v_012 = inv(PowerModelsProtection._A)*_v[5:7]
+@printf("secondary side transformer phase currents (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(i_line[1]), angle(i_line[1])*180/pi, abs(i_line[2]), angle(i_line[2])*180/pi, abs(i_line[3]), angle(i_line[3])*180/pi, abs(i_line[4]), angle(i_line[4])*180/pi)
+@printf("secondary side transformer phase voltages (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(_v[5]), angle(_v[5])*180/pi, abs(_v[6]), angle(_v[6])*180/pi, abs(_v[7]), angle(_v[7])*180/pi, abs(_v[8]), angle(_v[8])*180/pi)
+@printf("secondary side transformer sequence currents (0, 1 ,2) = %f @ %F, %f @ %f, %f @ %F\n", abs(i_012[1]), angle(i_012[1])*180/pi, abs(i_012[2]), angle(i_012[2])*180/pi,abs(i_012[3]), angle(i_012[3])*180/pi)
+@printf("secondary side transformer sequence voltages (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(v_012[1]), angle(v_012[1])*180/pi, abs(v_012[2]), angle(v_012[2])*180/pi,abs(v_012[3]), angle(v_012[3])*180/pi)
+
+
+
+"""
+    lg- fault at load bus
+"""
+
+gp = 1 / phase_resistance
+gf = 1 / ground_resistance
+Gf_lg = zeros(Real, 1, 1)
+
+Gf_lg[1,1] += gf
+
+
+                                                           
+bus = model.data["bus"]["3"]
+indx = (1,2)
+y = deepcopy(model.y)
+for i_indx in 1:1
+    for j_indx in 1:1
+        i = model.data["admittance_map"][(bus["bus_i"],bus["terminals"][i_indx])]
+        j = model.data["admittance_map"][(bus["bus_i"],bus["terminals"][j_indx])]
+        y[i,j] += Gf_lg[i_indx,j_indx]
+    end
+end
+
+v = compute_mc_pf(y, model)
+data = model.data
+gen = data["gen"]["1"]
+transformer = data["transformer"][gen["transformer_id"]]
+f_bus = data["bus"]["$(transformer["f_bus"])"]
+t_bus = data["bus"]["$(transformer["t_bus"])"]
+y = transformer["p_matrix"][5:8,1:8]
+_v = zeros(Complex{Float64}, 8, 1)
+indx = 1
+for (_j, j) in enumerate(f_bus["terminals"])
+    if haskey(data["admittance_map"], (f_bus["bus_i"], j))
+        _v[indx, 1] = v[data["admittance_map"][(f_bus["bus_i"], j)], 1]
+    else
+        _v[indx, 1] = 0.0
+    end
+    global indx += 1
+end
+for (_j, j) in enumerate(t_bus["terminals"])
+    if haskey(data["admittance_map"], (t_bus["bus_i"], j))
+        _v[indx, 1] = v[data["admittance_map"][(t_bus["bus_i"], j)], 1]
+    else
+        _v[indx, 1] = 0.0
+    end
+    global indx += 1
+end
+i_abc = (transformer["p_matrix"][1:4,1:8]*_v)
+i_line = (transformer["p_matrix"][5:8,1:8]*_v)
+i_012 = inv(PowerModelsProtection._A)*i_abc[1:3]
+v_012 = inv(PowerModelsProtection._A)*_v[1:3]
+println("_____lg-fault____")
+@printf("inverter phase currents (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(i_abc[1]), angle(i_abc[1])*180/pi, abs(i_abc[2]), angle(i_abc[2])*180/pi,abs(i_abc[3]), angle(i_abc[3])*180/pi)
+@printf("inverter phase voltages (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(_v[1]), angle(_v[1])*180/pi, abs(_v[2]), angle(_v[2])*180/pi, abs(_v[3]), angle(_v[3])*180/pi, abs(_v[4]), angle(_v[4])*180/pi)
+@printf("inverter sequence currents (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(i_012[1]), angle(i_012[1])*180/pi, abs(i_012[2]), angle(i_012[2])*180/pi,abs(i_012[3]), angle(i_012[3])*180/pi)
+@printf("inverter sequence voltages (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(v_012[1]), angle(v_012[1])*180/pi, abs(v_012[2]), angle(v_012[2])*180/pi,abs(v_012[3]), angle(v_012[3])*180/pi)
+
+i_012 = inv(PowerModelsProtection._A)*i_line[1:3]
+v_012 = inv(PowerModelsProtection._A)*_v[5:7]
+@printf("secondary side transformer phase currents (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(i_line[1]), angle(i_line[1])*180/pi, abs(i_line[2]), angle(i_line[2])*180/pi, abs(i_line[3]), angle(i_line[3])*180/pi, abs(i_line[4]), angle(i_line[4])*180/pi)
+@printf("secondary side transformer phase voltages (a, b, c) = %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f, %.2f @ %.2f\n", abs(_v[5]), angle(_v[5])*180/pi, abs(_v[6]), angle(_v[6])*180/pi, abs(_v[7]), angle(_v[7])*180/pi, abs(_v[8]), angle(_v[8])*180/pi)
+@printf("secondary side transformer sequence currents (0, 1 ,2) = %f @ %F, %f @ %f, %f @ %F\n", abs(i_012[1]), angle(i_012[1])*180/pi, abs(i_012[2]), angle(i_012[2])*180/pi,abs(i_012[3]), angle(i_012[3])*180/pi)
+@printf("secondary side transformer sequence voltages (0, 1, 2) = %f @ %F, %f @ %f, %f @ %F\n", abs(v_012[1]), angle(v_012[1])*180/pi, abs(v_012[2]), angle(v_012[2])*180/pi,abs(v_012[3]), angle(v_012[3])*180/pi)
\ No newline at end of file