From 432a796498eb47abf950b4820f8fb2cb1b305477 Mon Sep 17 00:00:00 2001
From: Changgang Li <lichgang@sdu.edu.cn>
Date: Mon, 22 Jul 2019 11:32:25 +0800
Subject: [PATCH] v0.3.23.2. (1) Fix bug in hvdc_model: the loss of inverter
 side should be deleted from the DC power. (2) Update the loss calcuation
 logic of HVDC device. (3) Add GENROU and GENSAL model checking function

---
 code/STEPS.depend                             | 12 ++---
 code/source/device/hvdc.cpp                   | 34 ++++++++++--
 code/source/model/hvdc_model/hvdc_model.cpp   |  9 +++-
 .../sg_models/sync_generator_model/genrou.cpp | 31 ++++++++++-
 .../sg_models/sync_generator_model/gensal.cpp | 21 +++++++-
 .../dynamic_simulator/dynamic_simulator.cpp   | 53 ++++++++++++++++++-
 6 files changed, 146 insertions(+), 14 deletions(-)

diff --git a/code/STEPS.depend b/code/STEPS.depend
index fa9af1fb..e32ff13f 100644
--- a/code/STEPS.depend
+++ b/code/STEPS.depend
@@ -2386,7 +2386,7 @@
 	"header/power_system_database.h"
 	"header/STEPS.h"
 
-1563500065 source:c:\users\lichg\steps\code\source\device\hvdc.cpp
+1563676767 source:c:\users\lichg\steps\code\source\device\hvdc.cpp
 	"header/device/hvdc.h"
 	"header/basic/utility.h"
 	"header/basic/constants.h"
@@ -2838,7 +2838,7 @@
 	"header/model/hvdc_model/hvdc_model_test.h"
 	"header/STEPS.h"
 
-1559112961 source:c:\users\lichg\steps\code\source\model\hvdc_model\hvdc_model.cpp
+1563700825 source:c:\users\lichg\steps\code\source\model\hvdc_model\hvdc_model.cpp
 	"header/model/hvdc_model/hvdc_model.h"
 	"header/steps_namespace.h"
 	"header/basic/utility.h"
@@ -3806,7 +3806,7 @@
 	"header/model/sg_models/sync_generator_model/gencls.h"
 	"header/power_system_database.h"
 
-1561097476 source:c:\users\lichg\steps\code\source\model\sg_models\sync_generator_model\genrou.cpp
+1563699881 source:c:\users\lichg\steps\code\source\model\sg_models\sync_generator_model\genrou.cpp
 	"header/model/sg_models/sync_generator_model/genrou.h"
 	"header/basic/utility.h"
 	"header/STEPS.h"
@@ -3834,7 +3834,7 @@
 	"header/model/sg_models/sync_generator_model/genrou.h"
 	"header/STEPS.h"
 
-1561097475 source:c:\users\lichg\steps\code\source\model\sg_models\sync_generator_model\gensal.cpp
+1563699929 source:c:\users\lichg\steps\code\source\model\sg_models\sync_generator_model\gensal.cpp
 	"header/model/sg_models/sync_generator_model/gensal.h"
 	"header/basic/utility.h"
 	"header/STEPS.h"
@@ -4798,7 +4798,7 @@
 	"header/toolkit/cct_searcher/cct_searcher.h"
 	"header/STEPS.h"
 
-1563500568 source:c:\users\lichg\steps\code\source\toolkit\dynamic_simulator\dynamic_simulator.cpp
+1563700780 source:c:\users\lichg\steps\code\source\toolkit\dynamic_simulator\dynamic_simulator.cpp
 	"header/toolkit/dynamic_simulator/dynamic_simulator.h"
 	"header/basic/utility.h"
 	"header/steps_namespace.h"
@@ -4841,7 +4841,7 @@
 	"header/toolkit/dynamic_simulator/dynamic_simulator.h"
 	<fstream>
 
-1563241935 source:c:\users\lichg\steps\code\source\toolkit\powerflow_solver\powerflow_solver.cpp
+1563678595 source:c:\users\lichg\steps\code\source\toolkit\powerflow_solver\powerflow_solver.cpp
 	"header/toolkit/powerflow_solver/powerflow_solver.h"
 	"header/basic/utility.h"
 	"header/steps_namespace.h"
diff --git a/code/source/device/hvdc.cpp b/code/source/device/hvdc.cpp
index 75f547e7..3d240a97 100644
--- a/code/source/device/hvdc.cpp
+++ b/code/source/device/hvdc.cpp
@@ -1182,6 +1182,12 @@ void HVDC::run(DYNAMIC_MODE mode)
         {
             case INITIALIZE_MODE:
             {
+                if(toolkit.is_detailed_log_enabled())
+                {
+                    osstream<<"HVDC data before dynamic initialization";
+                    toolkit.show_information_with_leading_time_stamp(osstream);
+                    show_solved_hvdc_steady_state();
+                }
                 HVDC_MODEL* hvdc_model = get_hvdc_model();
                 if(hvdc_model!=NULL)
                     hvdc_model->initialize();
@@ -1385,10 +1391,20 @@ void HVDC::show_solved_hvdc_steady_state() const
     osstream<<"Rectifier AC current = "<<get_converter_ac_current_in_kA(converter)<<" kA.";
     toolkit.show_information_with_leading_time_stamp(osstream);
 
+    double P = get_converter_ac_active_power_in_MW(converter);
+    double Q = get_converter_ac_reactive_power_in_MVar(converter);
+    complex<double> S(P,Q);
+    complex<double> V = psdb.get_bus_complex_voltage_in_pu(get_converter_bus(converter));
+    S /=  psdb.get_system_base_power_in_MVA();
+    complex<double> I = conj(S/V);
+    I *= (psdb.get_system_base_power_in_MVA()/(sqrt(3.0)*psdb.get_bus_base_voltage_in_kV(get_converter_bus(converter))));
+
     osstream<<"Pdc = "<<get_converter_dc_power_in_MW(converter)<<" MW, "
             <<"Pac = "<<get_converter_ac_active_power_in_MW(converter)<<" MW, "
             <<"Qac = "<<get_converter_ac_reactive_power_in_MVar(converter)<<" MVar, "
-            <<"Eta = "<<get_converter_ac_power_factor(converter);
+            <<"Eta = "<<get_converter_ac_power_factor(converter)<<", "
+            <<"Iac = "<<get_converter_ac_current_in_kA(converter)<<"kA, or "
+            <<I<<"KA and "<<abs(I)<<"kA";
     toolkit.show_information_with_leading_time_stamp(osstream);
 
     converter = INVERTER;
@@ -1438,10 +1454,20 @@ void HVDC::show_solved_hvdc_steady_state() const
     osstream<<"Inverter AC current = "<<get_converter_ac_current_in_kA(converter)<<" kA.";
     toolkit.show_information_with_leading_time_stamp(osstream);
 
+    P = get_converter_ac_active_power_in_MW(converter);
+    Q = get_converter_ac_reactive_power_in_MVar(converter);
+    S = complex<double>(P,Q);
+    V = psdb.get_bus_complex_voltage_in_pu(get_converter_bus(converter));
+    S /=  psdb.get_system_base_power_in_MVA();
+    I = conj(S/V);
+    I *= (psdb.get_system_base_power_in_MVA()/(sqrt(3.0)*psdb.get_bus_base_voltage_in_kV(get_converter_bus(converter))));
+
     osstream<<"Pdc = "<<get_converter_dc_power_in_MW(converter)<<" MW, "
             <<"Pac = "<<get_converter_ac_active_power_in_MW(converter)<<" MW, "
             <<"Qac = "<<get_converter_ac_reactive_power_in_MVar(converter)<<" MVar, "
-            <<"Eta = "<<get_converter_ac_power_factor(converter);
+            <<"Eta = "<<get_converter_ac_power_factor(converter)<<", "
+            <<"Iac = "<<get_converter_ac_current_in_kA(converter)<<"kA, or "
+            <<I<<"KA and "<<abs(I)<<"kA";
     toolkit.show_information_with_leading_time_stamp(osstream);
 }
 
@@ -1548,10 +1574,12 @@ double HVDC::get_converter_ac_active_power_in_MW(HVDC_CONVERTER_SIDE converter)
     return S*pf;*/
 
     double Iac = get_converter_ac_current_in_kA(converter);
+    double Idc = get_converter_dc_current_in_kA(converter);
 
     size_t N = get_converter_number_of_bridge(converter);
     double R = get_converter_transformer_impedance_in_ohm(converter).real();
-    double loss = sqrt(3.0)*Iac*Iac*R/N;
+    //double loss = sqrt(3.0)*Iac*Iac*R/N;
+    double loss = 2.0*Idc*Idc*R*N;
 
     switch(converter)
     {
diff --git a/code/source/model/hvdc_model/hvdc_model.cpp b/code/source/model/hvdc_model/hvdc_model.cpp
index e9766274..08290f2f 100644
--- a/code/source/model/hvdc_model/hvdc_model.cpp
+++ b/code/source/model/hvdc_model/hvdc_model.cpp
@@ -818,7 +818,7 @@ void HVDC_MODEL::solve_hvdc_model_without_line_dynamics(double Iset_kA, double V
     double TIME = toolkit.get_dynamic_simulation_time_in_s();
 
     //osstream<<"solving "<<get_device_name()<<" with I command = "<<Iset_kA<<"kA, V command = "<<Vset_kV;
-    //show_information_with_leading_time_stamp(osstream);
+    //toolkit.show_information_with_leading_time_stamp(osstream);
 
     POWER_SYSTEM_DATABASE& psdb = toolkit.get_power_system_database();
 
@@ -863,6 +863,8 @@ void HVDC_MODEL::solve_hvdc_model_without_line_dynamics(double Iset_kA, double V
 
     double vac_r = psdb.get_bus_voltage_in_kV(bus_r);
     double vac_i = psdb.get_bus_voltage_in_kV(bus_i);
+    //osstream<<"rectifier and inverter side AC voltages are: "<<psdb.get_bus_voltage_in_pu(bus_r)<<", "<<psdb.get_bus_voltage_in_pu(bus_i);
+    //toolkit.show_information_with_leading_time_stamp(osstream);
 
     double eac_r = vac_r/tap_r*ebase_converter_r/ebase_grid_r;
     double eac_i = vac_i/tap_i*ebase_converter_i/ebase_grid_i;
@@ -1415,7 +1417,10 @@ complex<double> HVDC_MODEL::get_converter_ac_complex_power_in_MVA(HVDC_CONVERTER
     double Idc = get_converter_dc_current_in_kA(converter);
     double P = get_converter_dc_power_in_MW(converter);
     double Ploss = N*(2.0*Idc*Idc*Rc); // N for bridges, 2 for two phases
-    P += Ploss;
+    if(converter==RECTIFIER)
+        P += Ploss;
+    else
+        P -= Ploss;
     double phi = deg2rad(get_converter_ac_power_factor_angle_in_deg(converter));
     if(phi==0.0)
         return P;
diff --git a/code/source/model/sg_models/sync_generator_model/genrou.cpp b/code/source/model/sg_models/sync_generator_model/genrou.cpp
index 426a6b40..58eaccfc 100644
--- a/code/source/model/sg_models/sync_generator_model/genrou.cpp
+++ b/code/source/model/sg_models/sync_generator_model/genrou.cpp
@@ -484,7 +484,36 @@ double GENROU::get_terminal_current_in_pu_based_on_sbase()
 
 void GENROU::check()
 {
-    ;
+    ostringstream osstream;
+    STEPS& toolkit = get_toolkit(__PRETTY_FUNCTION__);
+    double xdp = get_Xdp();
+    double xqp = get_Xqp();
+    double xpp = get_Xpp();
+    double xl = get_Xl();
+    osstream<<"Error is detected at "<<get_model_name()<<" model of "<<get_device_name()<<".\n";
+    bool error_found = false;
+    if(xdp==xpp)
+    {
+        osstream<<"Xd'=X\" was detected: Xd'="<<xdp<<", X\"="<<xpp<<"\n";
+        error_found = true;
+    }
+    if(xqp==xpp)
+    {
+        osstream<<"Xq'=X\" was detected: Xq'="<<xqp<<", X\"="<<xpp<<"\n";
+        error_found = true;
+    }
+    if(xdp==xl)
+    {
+        osstream<<"Xd'=Xl was detected: Xd'="<<xdp<<", Xl="<<xl<<"\n";
+        error_found = true;
+    }
+    if(xqp==xl)
+    {
+        osstream<<"Xq'=Xl was detected: Xq'="<<xqp<<", Xl="<<xl<<"\n";
+        error_found = true;
+    }
+    if(error_found)
+        toolkit.show_information_with_leading_time_stamp(osstream);
 }
 
 void GENROU::report()
diff --git a/code/source/model/sg_models/sync_generator_model/gensal.cpp b/code/source/model/sg_models/sync_generator_model/gensal.cpp
index d3143f03..fa81370c 100644
--- a/code/source/model/sg_models/sync_generator_model/gensal.cpp
+++ b/code/source/model/sg_models/sync_generator_model/gensal.cpp
@@ -417,7 +417,26 @@ double GENSAL::get_terminal_current_in_pu_based_on_sbase()
 
 void GENSAL::check()
 {
-    ;
+    ostringstream osstream;
+    STEPS& toolkit = get_toolkit(__PRETTY_FUNCTION__);
+    double xdp = get_Xdp();
+    double xqp = get_Xqp();
+    double xpp = get_Xpp();
+    double xl = get_Xl();
+    osstream<<"Error is detected at "<<get_model_name()<<" model of "<<get_device_name()<<".\n";
+    bool error_found = false;
+    if(xdp==xpp)
+    {
+        osstream<<"Xd'=X\" was detected: Xd'="<<xdp<<", X\"="<<xpp<<"\n";
+        error_found = true;
+    }
+    if(xdp==xl)
+    {
+        osstream<<"Xd'=Xl was detected: Xd'="<<xdp<<", Xl="<<xl<<"\n";
+        error_found = true;
+    }
+    if(error_found)
+        toolkit.show_information_with_leading_time_stamp(osstream);
 }
 
 void GENSAL::report()
diff --git a/code/source/toolkit/dynamic_simulator/dynamic_simulator.cpp b/code/source/toolkit/dynamic_simulator/dynamic_simulator.cpp
index bdccfc53..827899f2 100644
--- a/code/source/toolkit/dynamic_simulator/dynamic_simulator.cpp
+++ b/code/source/toolkit/dynamic_simulator/dynamic_simulator.cpp
@@ -16,7 +16,7 @@ DYNAMICS_SIMULATOR::DYNAMICS_SIMULATOR()
 {
     DELT = 0.01;
     clear();
-    check_NAN = false;
+    check_NAN = true;
 }
 
 DYNAMICS_SIMULATOR::~DYNAMICS_SIMULATOR()
@@ -2104,6 +2104,26 @@ void DYNAMICS_SIMULATOR::get_bus_current_mismatch()
     //#pragma omp parallel for
     for(size_t i = 0; i<n; ++i)
         I_mismatch[i] = -I_mismatch[i];
+    if(check_NAN)
+    {
+        for(size_t i = 0; i<n; ++i)
+        {
+            if(isnan(I_mismatch[i].real()) or isnan(I_mismatch[i].imag()))
+            {
+                osstream<<"warning. NAN is detected when getting bus current mismatch after get_bus_currnet_into_network() when calling DYNAMICS_SIMULATOR::"<<__FUNCTION__<<"():"<<endl;
+                for(size_t j = 0; j<n; ++j)
+                {
+                    if(isnan(I_mismatch[j].real()) or isnan(I_mismatch[j].imag()))
+                    {
+                        size_t ibus = net.get_physical_bus_number_of_internal_bus(j);
+                        osstream<<"Physical bus: "<<ibus<<", internal bus: "<<j<<", "<<I_mismatch[i].real()<<","<<I_mismatch[i].imag()<<endl;
+                    }
+                }
+                toolkit.show_information_with_leading_time_stamp(osstream);
+                break;
+            }
+        }
+    }
 
     add_generators_to_bus_current_mismatch();
     if(check_NAN)
@@ -2231,6 +2251,20 @@ void DYNAMICS_SIMULATOR::get_bus_current_mismatch()
     }
     osstream<<"max mismatch @ bus "<<busmax<<", "<<maxmismatch<<endl;
     toolkit.show_information_with_leading_time_stamp(osstream);
+
+    for(size_t i=0; i<n; ++i)
+    {
+        if(abs(I_mismatch[i])>0.0001)
+        {
+            osstream<<"mismatch @ bus "<<net.get_physical_bus_number_of_internal_bus(i)<<", "<<I_mismatch[i]<<">0.0001"<<endl;
+            toolkit.show_information_with_leading_time_stamp(osstream);
+        }
+        if(isnan(I_mismatch[i].real()) or isnan(I_mismatch[i].imag()))
+        {
+            osstream<<"NAN is detected in I_mismatch @ bus "<<net.get_physical_bus_number_of_internal_bus(i)<<", "<<I_mismatch[i]<<endl;
+            toolkit.show_information_with_leading_time_stamp(osstream);
+        }
+    }
     */
 }
 
@@ -2402,6 +2436,7 @@ void DYNAMICS_SIMULATOR::add_loads_to_bus_current_mismatch()
 }
 void DYNAMICS_SIMULATOR::add_hvdcs_to_bus_current_mismatch()
 {
+    ostringstream osstream;
     STEPS& toolkit = get_toolkit(__PRETTY_FUNCTION__);
     POWER_SYSTEM_DATABASE& psdb = toolkit.get_power_system_database();
     NETWORK_MATRIX& network_matrix = get_network_matrix();
@@ -2427,6 +2462,14 @@ void DYNAMICS_SIMULATOR::add_hvdcs_to_bus_current_mismatch()
 
             I_mismatch[internal_bus] -= I;
 
+            if(toolkit.is_detailed_log_enabled())
+            {
+                I *= (psdb.get_system_base_power_in_MVA()/(sqrt(3.0)*psdb.get_bus_base_voltage_in_kV(physical_bus)));
+                osstream<<"Current at rectifier side of "<<hvdc->get_device_name()<<": "<<I<<"kA or "<<abs(I)<<"kA"<<endl
+                        <<"Complex power = "<<sqrt(3.0)*psdb.get_bus_complex_voltage_in_kV(physical_bus)*conj(I)<<" MVA";
+                toolkit.show_information_with_leading_time_stamp(osstream);
+            }
+
             physical_bus = hvdcs[i]->get_converter_bus(INVERTER);
 
             internal_bus = network_matrix.get_internal_bus_number_of_physical_bus(physical_bus);
@@ -2434,6 +2477,14 @@ void DYNAMICS_SIMULATOR::add_hvdcs_to_bus_current_mismatch()
             I = hvdcs[i]->get_converter_dynamic_current_in_pu_based_on_system_base_power(INVERTER);
 
             I_mismatch[internal_bus] -= I;
+
+            if(toolkit.is_detailed_log_enabled())
+            {
+                I *= (psdb.get_system_base_power_in_MVA()/(sqrt(3.0)*psdb.get_bus_base_voltage_in_kV(physical_bus)));
+                osstream<<"Current at inverter side of "<<hvdc->get_device_name()<<": "<<I<<"kA or "<<abs(I)<<"kA"<<endl
+                        <<"Complex power = "<<sqrt(3.0)*psdb.get_bus_complex_voltage_in_kV(physical_bus)*conj(I)<<" MVA";
+                toolkit.show_information_with_leading_time_stamp(osstream);
+            }
         }
     }
 }