diff --git a/EMTP-RV/ADN EPEC2017/DFIG_mask.dwj b/EMTP-RV/ADN EPEC2017/DFIG_mask.dwj new file mode 100644 index 0000000..5972f8d --- /dev/null +++ b/EMTP-RV/ADN EPEC2017/DFIG_mask.dwj @@ -0,0 +1,201 @@ +// Last change: HANI SAAD - JEAN MAHSEREDJIAN 04/11/2013 03:36:19 AM +//*This script is used for the DFIG wind generator model. +// This scripts allows to modify parameters that are not accessible on the main mask. + +with (Math) { + + //////////////////////////////////////////////////_BASIC CALCULATIONS ////////////////////////////////////////////////// + //*Fundamental Period + Period = 1 / Freq; + period = Period; + w = 2 * PI * Freq; + //*Total Wind park S and P + Pgen_W = Pgen_MW * 1e6; // Rated active power of wind generator (W) + Sgen_VA = Pgen_W * 1.11; // Rated apparent power of wind generator (VA). Note that 11% of reserve power is added. + Spark_VA = Ngen * Sgen_VA; // Total aggregated wind park S (VA) + Spark_MVA = Ngen * Sgen_VA / 1e6; // Total aggregated wind park S (MVA), ASM machine uses MVA selection + Ppark_W = Ngen * Pgen_W; // Total aggregated wind park P (W) (used for the load-flow device) + Qpark_Var = Qref_pu * Spark_VA; // Total aggregated wind park Q (Var) (used for the load-flow device when Q or Q+Vac control is selected) + Vdc_V = Vdc_kV * 1e3; // DC voltage (V) + //*base calculations + Vbase_VRMSLL = Vgen_kVRMSLL * 1e3; // Rated generator voltage (VRMSLL) + Vbase_V = Vbase_VRMSLL * sqrt(2 / 3); // Peak generator voltage (V) + Ibase_A = Spark_VA / Vbase_VRMSLL * sqrt(2 / 3); // Peak generator current (A) + Zbase_Ohm = (Vgen_kVRMSLL * Vgen_kVRMSLL) / Spark_MVA; // Base impedance (Ohm) + + //////////////////////////////////////////////////_ASYNCHRONOUS GENERATOR & AC/DC/AC CONVERTER PARAMETERS ////////////////////////////////////////////////// + //*ASM parameters (see ASM_DFIG/ASM1) + npoles = 3; // KTH, original value is 6 Number of poles of the asynchronous generator + Rs_pu = 0.00706 // KTH, original value is 0.007 in pu + Lls_pu = 0.171; // in pu + Lmd_pu = 2.9 // in pu + Lmq_pu = 2.9 // in pu + Rr_pu = 0.005 // in pu + Llr_pu = 0.156 // in pu + Hgen_pu = 5.06 // KTH=5.04 in sec + D_pu = 0.0 // in pu + Llr_plus_Lmd_pu = Llr_pu + Lmd_pu; // in pu (used in DFIG_CONTROL/ROTOR side CONTROL) + //*Choke impedance (see ASM_DFIG/Choke) + Rchoke_pu = 0.003; // in pu + Rchoke_Ohm = 1.8034e-5; // KTH, original value is "Rchoke_pu * Zbase_Ohm " in Ohm + Lchoke_pu = 0.3 // in pu + Lchoke_H = 5.7404e-6; // KTH, original value is "Lchoke_pu * Zbase_Ohm / (2 * PI * Freq)" in H + //*AC/DC/AC Converter parameters (see ASM_DFIG/BtoB_converter) + E_Cdc = 5; // in KJ/MVA + Cdc_F = 0.01; // KTH=0.01, original value is "E_Cdc * 1e-3 / (1 / 2 * pow(Vdc_V, 2) / Sgen_VA); " DC capacitor (F) + CdcPark_F = Cdc_F * Ngen; // Aggregated DC capacitor (F) + CarrierSignal_Freq = 3e3; // PWM Carrier signal frequency (Hz) (used only for Detailed Converter Model) + CarrierSignal_ratio = CarrierSignal_Freq / Freq; // PWM carrier signal ratio (used only for Detailed Converter Model) + Rsnubber_Ohm = 3000; // IGBT/diode resistance snubber (Ohm) (used only for Detailed Converter Model) + Csnubber_F = 1e-6; // IGBT/diode capacitor snubber (F) (used only for Detailed Converter Model) + Rchopper_Ohm = 0.2; // Chopper resistance in Ohm + RchopperPark_Ohm = Rchopper_Ohm / Ngen; // Aggregation of the chopper resistance + //*High pass filter parameters (see ASM_DFIG/filter_stator_ASM ) + Qfilter_Var = 0.473e6 // KTH setting, original value is "150e3 * Ngen" Reactive power generated by each filter + n1 = CarrierSignal_ratio // First tuning harmonic + n2 = CarrierSignal_ratio * 2 // Second tuning harmonic + QF = 1e3 // Quality factor + Cfilter1_F = Qfilter_Var / (Vbase_VRMSLL * Vbase_VRMSLL * 2 * PI * Freq); + Lfilter1_H = 1e6; // KTH, original value is " 1 / (Cfilter1_F * pow((n1 * Freq * 2 * PI), 2));" + Rfilter1_Ohm = 1.7 // KTH, original value is "QF * 2 * PI * n1 * Freq * Lfilter1_H;" + Cfilter2_F = Cfilter1_F; + Lfilter2_H = 1 / (Cfilter2_F * pow((n2 * Freq * 2 * PI), 2)); + Rfilter2_Ohm = QF * 2 * PI * n2 * Freq * Lfilter2_H; + + //////////////////////////////////////////////////_WIND TURBINE PARAMETERS ////////////////////////////////////////////////// + //* Blade and Gearbox parameters (see WIND_TURBINE/Wind_Power) + blade_radius = ceil((sqrt(pow(35.25, 2) / 1.5e6 * Pgen_W)) * 100) / 100; // Blade radius (m) is in function of the rated power + gearbox_ratio = 72 * (Freq / 60); // Gearbox ratio. Note that the ratio is adjusted depending on the frequency used in order to preserve the same tip-speed value (Lambda). + ro = 1.225; // air density (kg/m^3) + K_speed = 2 * PI * Freq / (npoles / 2) * blade_radius / gearbox_ratio; // speed coefficient defined as: Lambda = K_speed / wind_speed * w_rotor; + K_air = (1 / 2) * PI * blade_radius * blade_radius * ro; // air coefficient defined as: Pair = K_air * Cp * (wind_speed)^3; + Pitch_max_deg = 45; // KTH, original value is 20 Maximum pitch angle in degree + Pitch_min_deg = 0; // KTH, original value is -1.5 Minimum pitch angle in degree + + //*Shaft and train drive parameters (see WIND_TURBINE/Two_mass_model) + Ksh = 1.2; //Speed deviation damping + Dm = 1.5; //Mutual damping + Ht = 4 // Inertia constant (s) + + //////////////////////////////////////////////////_TRANSFORMER PARAMETERS ////////////////////////////////////////////////// + //* Wind generator transformer parameters (see Transformer) + Rtransfo_sec_Ohm = (Rtransfo_pu * 0.1 * Vgrid_kVRMSLL * Vgrid_kVRMSLL / Spark_MVA); // Secondary winding resistance, Ohm + Ltransfo_sec_H = (Xtransfo_pu * 0.1 * Vgrid_kVRMSLL * Vgrid_kVRMSLL / Spark_MVA) / (2 * PI * Freq); // Secondary winding inductance, H + Rtransfo_prim_Ohm = (Rtransfo_pu * 0.9 * Vgen_kVRMSLL * Vgen_kVRMSLL / Spark_MVA); // Primary winding resistance, Ohm + Ltransfo_prim_H = (Xtransfo_pu * 0.9 * Vgen_kVRMSLL * Vgen_kVRMSLL / Spark_MVA) / (2 * PI * Freq); // Primary winding inductance, H + Transfo_ratio = Vgen_kVRMSLL / Vgrid_kVRMSLL / sqrt(3); // Turns ratio + + //////////////////////////////////////////////////_CONTROL SYSTEM PARAMETERS ////////////////////////////////////////////////// + //* Pitch control (see CONTROL/Pitch_Control ) + // NB: Ptich Control Gains are tuned by trial and error + PitchCtrl_ki = 20; // Integral Gain of the Pitch control + PitchCtrl_kp = 500; // KTH, original value is 50 Proportional Gain of the Pitch control + PitchCtrl_Td_d = 10; // Derivative Gain of the Pitch control + PitchCtrl_tau_d = 5e-3; // derivative filter in s + w_rotor_max_pu = 1.25; // maximum rotor speed in pu + Pitch_rate = 2 // KTH, original value is 10 Slop rate (see Pitch_Control/Rate_Limiter) + + //*Rotor control (see CONTROL/Rotor_Control ) + //+Rotor Inner current control + eps = 0.7; // Damping ratio + w_i = 3 / (eps * TimeConstant_i_rotor); // Impulse of the closed loop control + RotorCtrl_ki = 8 // KTH, 0.2 * ((Llr_plus_Lmd_pu) / w) * w_i * w_i; Integral Gain of the Inner Current control + RotorCtrl_kp = 0.3; // KTH, original value is 0.2 * 2 * eps * w_i * ((Llr_plus_Lmd_pu) / w); Proportional gain of the Inner Current control + //+Q control and Q+Vac control + Qctrl_kp = 0.05; //KTH, original value is 0 Proportional Gain of Q control (Only Integral part of the PI control is used) + Qctrl_ki = 5; // KTH, original value is "3 / TimeConstant_Q" Integral Gain of Q control + //+Vac control + X_tot = X_PCC_pu + Xtransfo_pu; // Total impedance in pu at wind generator connection point + VacCtrl_ki = 3 / (X_tot * TimeConstant_Vac); // Integral Gain of Vac control + VacCtrl_kp = 1; // Proportional Gain of Vac control (Only Integral part of the PI control is used) + //+P control + Pctrl_kp = 1; // Proportional Gain of P control (Only Integral part of the PI control is used) + Pctrl_ki = 3 / TimeConstant_P; // Integral Gain of P control + //+Maximum Power Point Tracking (MPPT) + Lambda_optimal = 6.635; // Optimal Lambda in order to have the rated active power for a wind speed = 11.24 m/s + Cp_max = 0.46; // Related Cp_max in order to have the rated active power for a wind speed = 11.24 m/s + K_optimal_pu = K_air * Cp_max * pow(K_speed / Lambda_optimal, 3) / Sgen_VA; // Constant coefficient defined as: P = K_optimal_pu * w^3 + //+ Idq limit (see CONTROL/Rotor_Control/Idq_ref_limiter1 ) + Q_priority = 1; // 1 = priority is given for reactive power (should be set for FRT capability). 0 = priority is given to active power + I_lim_Rotor_pu = 1.1; // Maximum rotor Idq limit + Id_lim_Rotor_pu = 1; // Maximum rotor Id limit, if Q_priority = 1 + Iq_lim_Rotor_pu = 1.1; // Maximum rotor Iq limit, if Q_priority = 0 + + //*Grid control (see CONTROL/Grid_Control ) + //+Grid Inner current control + eps = 0.7; // Damping ratio + w_i = 3 / (eps * TimeConstant_i_grid); // Impulse of the closed loop control + GridCtrl_ki = 500; // KTH, original value is "(Lchoke_pu / w) * w_i * w_i" Integral Gain of the Inner Current control + GridCtrl_kp = 2.5; // KTH, original value is "2 * eps * w_i * (Lchoke_pu / w)" Proportional gain of the Inner Current control + //+Vdc control + E_Cdc_J = 0.5 * CdcPark_F * Vdc_V * Vdc_V; // Energy stored in the DC capacitor (Joule) + H_Cdc = E_Cdc_J / (0.3 * Spark_VA); // Static moment of inertia (S). Note: Power converter is rated for 30% of the nominal power + eps = 1; // Damping ratio + w_Vdc = 3 / (eps * TimeConstant_Vdc); // Time response for 5% (sec) + VdcCtrl_ki = H_Cdc * w_Vdc * w_Vdc; // KTH=0.05, Proportional Gain of Vdc control + VdcCtrl_kp = 2 * eps * w_Vdc * H_Cdc; // KTH=2e-3 Integral Gain of Vdc control + + //*PLL (see CONTROL/Computing_Variables/PLL) + PLL_lim_low = 12 + PLL_lim_up = 12 + PLL_max_delay = 0.0223 + eps = 1; // damping ratio + w_i = 3 / (eps * TimeConstant_PLL); + PLL_ki = 1 * w_i * w_i; // Int. Gain integral (rad/V) + PLL_kp = 2 * eps * w_i * 1; // Prop. Gain (rad/s/V) + + //*Protection system (see CONTROL/Protection ) + //+Deep Voltage Sag (DVS) protection + DVS_level = 0.001; // Threshold Voltage value in pu to activate the DVS protection + DVS_hysteresis = 0.1; // Threshold Voltage value in pu to deactivate the DVS protection + //+Crowbar protection + Rcrowbar_Ohm = 0.09; // Crowbar resistance (Ohm) + Lcrowbar_H = 4.0e-05; // Crowbar inductance (H) + RcrowbarPark_Ohm = Rcrowbar_Ohm / Ngen; // Crowbar equivalent resistance (Ohm) + LcrowbarPark_H = Lcrowbar_H / Ngen; // Crowbar equivalent inductance (H) + //+Converter Overcurrent limit + I_RotorSideConv_max = 4; // Threshold Current value in pu to activate RotorSide Overcurrent protection + I_GridSideConv_max = 4; // Threshold Current value in pu to activate GridSide Overcurrent protection + + //////////////////////////////////////////////////_INITIALIZATION ////////////////////////////////////////////////// + //*Initialization of the Protection system (see CONTROL/Protection) + init_Prot_delay = 0.3; // Delay time initialization of the Protection devices (s) (All protections, expect the Chopper, are deactivated for the specified time ) + + //*Initialization of the WIND TURBINE model (see WIND_TURBINE/Two_mass_model) + Lambda_ini_guess = Lambda_optimal; // Initial Lambda guess value + Cp_init_guess = Cp_max; // Initial Cp guess value + max_wind_speed = K_speed * w_rotor_max_pu / Lambda_ini_guess; // maximum wind speed before the Pitch control is activated + if (Mean_wind_speed > max_wind_speed) { // if the wind speed value specified by the user is higher than the maximum wind speed before Pitch control is activated + Mean_wind_speed_init = max_wind_speed; // than the wind speed value is fixed at his maximum value for initialization. + } + else { + Mean_wind_speed_init = Mean_wind_speed; + } + + init_w_turbine_pu = Mean_wind_speed_init * Lambda_ini_guess / K_speed * 0.989; // Initial Turbine speed in pu. The value is multiplied by a factor of losses. + init_Pturbine_pu = K_air * Cp_init_guess * pow(Mean_wind_speed_init, 3) / Sgen_VA; + init_Tturbine_pu = init_Pturbine_pu / init_w_turbine_pu; // Initial Turbine Torque in pu. + + //*Initialization of _ASYNCHRONOUS GENERATOR model (See ASM_DFIG/ASM1) + init_slip_pct = (1 - init_w_turbine_pu / 1) * 100; // Initial slip (%) + + /////////////////////////////////////////////////////////////_LOAD-FLOW CALCULATIONS ////////////////////////////////////////////////// + + //* Load-Flow of the wind generators (see LF_init1) + LF_V = Vgen_kVRMSLL * 1000 * sqrt(2 / 3); // Voltage for the LF device at the bus of the generator + + if (ControlMode == 1 || ControlMode == 3) { // if Q control or Q+Vac Control is selected, than PQ constraint is activated + LF_t1 = -1; // switch trick to active PQ constraint + LF_t2 = 1e15; // switch trick to active PQ constraint + LF_Ppark1_W = init_Pturbine_pu * Spark_VA * 0.967; // Active power in W of PQ constraint. The value is multiplied by a factor of losses. + LF_Ppark2_W = 0; // Active power in W of PV constraint. The value is set to zero + } else { // if Vac control is selected, than PV constraint is activated + LF_t1 = 1e15; // switch trick to active PV constraint + LF_t2 = -1; // switch trick to active PV constraint + LF_Ppark1_W = 0; // Active power in W of PQ constraint. The value is set to zero + LF_Ppark2_W = init_Pturbine_pu * Spark_VA * 0.967; // Active power in W of PV constraint. The value is multiplied by a factor of losses. + } + + +} + diff --git a/EMTP-RV/ADN EPEC2017/DFIG_mask2.dwj b/EMTP-RV/ADN EPEC2017/DFIG_mask2.dwj new file mode 100644 index 0000000..7709f84 --- /dev/null +++ b/EMTP-RV/ADN EPEC2017/DFIG_mask2.dwj @@ -0,0 +1,201 @@ +// Last change: HANI SAAD - JEAN MAHSEREDJIAN 04/11/2013 03:36:19 AM +//*This script is used for the DFIG wind generator model. +// This scripts allows to modify parameters that are not accessible on the main mask. + +with (Math) { + + //////////////////////////////////////////////////_BASIC CALCULATIONS ////////////////////////////////////////////////// + //*Fundamental Period + Period = 1 / Freq; + period = Period; + w = 2 * PI * Freq; + //*Total Wind park S and P + Pgen_W = Pgen_MW * 1e6; // Rated active power of wind generator (W) + Sgen_VA = Pgen_W * 1.11; // Rated apparent power of wind generator (VA). Note that 11% of reserve power is added. + Spark_VA = Ngen * Sgen_VA; // Total aggregated wind park S (VA) + Spark_MVA = Ngen * Sgen_VA / 1e6; // Total aggregated wind park S (MVA), ASM machine uses MVA selection + Ppark_W = Ngen * Pgen_W; // Total aggregated wind park P (W) (used for the load-flow device) + Qpark_Var = Qref_pu * Spark_VA; // Total aggregated wind park Q (Var) (used for the load-flow device when Q or Q+Vac control is selected) + Vdc_V = Vdc_kV * 1e3; // DC voltage (V) + //*base calculations + Vbase_VRMSLL = Vgen_kVRMSLL * 1e3; // Rated generator voltage (VRMSLL) + Vbase_V = Vbase_VRMSLL * sqrt(2 / 3); // Peak generator voltage (V) + Ibase_A = Spark_VA / Vbase_VRMSLL * sqrt(2 / 3); // Peak generator current (A) + Zbase_Ohm = (Vgen_kVRMSLL * Vgen_kVRMSLL) / Spark_MVA; // Base impedance (Ohm) + + //////////////////////////////////////////////////_ASYNCHRONOUS GENERATOR & AC/DC/AC CONVERTER PARAMETERS ////////////////////////////////////////////////// + //*ASM parameters (see ASM_DFIG/ASM1) + npoles = 6; // Number of poles of the asynchronous generator + Rs_pu = 0.007 // in pu + Lls_pu = 0.171; // in pu + Lmd_pu = 2.9 // in pu + Lmq_pu = 2.9 // in pu + Rr_pu = 0.005 // in pu + Llr_pu = 0.156 // in pu + Hgen_pu = 5.04 // in sec + D_pu = 0.0 // in pu + Llr_plus_Lmd_pu = Llr_pu + Lmd_pu; // in pu (used in DFIG_CONTROL/ROTOR side CONTROL) + //*Choke impedance (see ASM_DFIG/Choke) + Rchoke_pu = 0.003; // in pu + Rchoke_Ohm = Rchoke_pu * Zbase_Ohm; // in Ohm + Lchoke_pu = 0.3 // in pu + Lchoke_H = Lchoke_pu * Zbase_Ohm / (2 * PI * Freq); // in H + //*AC/DC/AC Converter parameters (see ASM_DFIG/BtoB_converter) + E_Cdc = 5; // in KJ/MVA + Cdc_F = E_Cdc * 1e-3 / (1 / 2 * pow(Vdc_V, 2) / Sgen_VA); // DC capacitor (F) + CdcPark_F = Cdc_F * Ngen; // Aggregated DC capacitor (F) + CarrierSignal_Freq = 3e3; // PWM Carrier signal frequency (Hz) (used only for Detailed Converter Model) + CarrierSignal_ratio = CarrierSignal_Freq / Freq; // PWM carrier signal ratio (used only for Detailed Converter Model) + Rsnubber_Ohm = 3000; // IGBT/diode resistance snubber (Ohm) (used only for Detailed Converter Model) + Csnubber_F = 1e-6; // IGBT/diode capacitor snubber (F) (used only for Detailed Converter Model) + Rchopper_Ohm = 0.2; // Chopper resistance in Ohm + RchopperPark_Ohm = Rchopper_Ohm / Ngen; // Aggregation of the chopper resistance + //*High pass filter parameters (see ASM_DFIG/filter_stator_ASM ) + Qfilter_Var = 150e3 * Ngen // Reactive power generated by each filter + n1 = CarrierSignal_ratio // First tuning harmonic + n2 = CarrierSignal_ratio * 2 // Second tuning harmonic + QF = 1e3 // Quality factor + Cfilter1_F = Qfilter_Var / (Vbase_VRMSLL * Vbase_VRMSLL * 2 * PI * Freq); + Lfilter1_H = 1 / (Cfilter1_F * pow((n1 * Freq * 2 * PI), 2)); + Rfilter1_Ohm = QF * 2 * PI * n1 * Freq * Lfilter1_H; + Cfilter2_F = Cfilter1_F; + Lfilter2_H = 1 / (Cfilter2_F * pow((n2 * Freq * 2 * PI), 2)); + Rfilter2_Ohm = QF * 2 * PI * n2 * Freq * Lfilter2_H; + + //////////////////////////////////////////////////_WIND TURBINE PARAMETERS ////////////////////////////////////////////////// + //* Blade and Gearbox parameters (see WIND_TURBINE/Wind_Power) + blade_radius = ceil((sqrt(pow(35.25, 2) / 1.5e6 * Pgen_W)) * 100) / 100; // Blade radius (m) is in function of the rated power + gearbox_ratio = 72 * (Freq / 60); // Gearbox ratio. Note that the ratio is adjusted depending on the frequency used in order to preserve the same tip-speed value (Lambda). + ro = 1.225; // air density (kg/m^3) + K_speed = 2 * PI * Freq / (npoles / 2) * blade_radius / gearbox_ratio; // speed coefficient defined as: Lambda = K_speed / wind_speed * w_rotor; + K_air = (1 / 2) * PI * blade_radius * blade_radius * ro; // air coefficient defined as: Pair = K_air * Cp * (wind_speed)^3; + Pitch_max_deg = 20; // Maximum pitch angle in degree + Pitch_min_deg = -1.5; // Minimum pitch angle in degree + + //*Shaft and train drive parameters (see WIND_TURBINE/Two_mass_model) + Ksh = 1.2; //Speed deviation damping + Dm = 1.5; //Mutual damping + Ht = 4.0 // Inertia constant (s) + + //////////////////////////////////////////////////_TRANSFORMER PARAMETERS ////////////////////////////////////////////////// + //* Wind generator transformer parameters (see Transformer) + Rtransfo_sec_Ohm = (Rtransfo_pu * 0.1 * Vgrid_kVRMSLL * Vgrid_kVRMSLL / Spark_MVA); // Secondary winding resistance, Ohm + Ltransfo_sec_H = (Xtransfo_pu * 0.1 * Vgrid_kVRMSLL * Vgrid_kVRMSLL / Spark_MVA) / (2 * PI * Freq); // Secondary winding inductance, H + Rtransfo_prim_Ohm = (Rtransfo_pu * 0.9 * Vgen_kVRMSLL * Vgen_kVRMSLL / Spark_MVA); // Primary winding resistance, Ohm + Ltransfo_prim_H = (Xtransfo_pu * 0.9 * Vgen_kVRMSLL * Vgen_kVRMSLL / Spark_MVA) / (2 * PI * Freq); // Primary winding inductance, H + Transfo_ratio = Vgen_kVRMSLL / Vgrid_kVRMSLL / sqrt(3); // Turns ratio + + //////////////////////////////////////////////////_CONTROL SYSTEM PARAMETERS ////////////////////////////////////////////////// + //* Pitch control (see CONTROL/Pitch_Control ) + // NB: Ptich Control Gains are tuned by trial and error + PitchCtrl_ki = 20; // Integral Gain of the Pitch control + PitchCtrl_kp = 50; // Proportional Gain of the Pitch control + PitchCtrl_Td_d = 10; // Derivative Gain of the Pitch control + PitchCtrl_tau_d = 5e-3; // derivative filter in s + w_rotor_max_pu = 1.25; // maximum rotor speed in pu + Pitch_rate = 10 // Slop rate (see Pitch_Control/Rate_Limiter) + + //*Rotor control (see CONTROL/Rotor_Control ) + //+Rotor Inner current control + eps = 0.7; // Damping ratio + w_i = 3 / (eps * TimeConstant_i_rotor); // Impulse of the closed loop control + RotorCtrl_ki = 0.2 * ((Llr_plus_Lmd_pu) / w) * w_i * w_i; // Integral Gain of the Inner Current control + RotorCtrl_kp = 0.2 * 2 * eps * w_i * ((Llr_plus_Lmd_pu) / w); // Proportional gain of the Inner Current control + //+Q control and Q+Vac control + Qctrl_kp = 0; // Proportional Gain of Q control (Only Integral part of the PI control is used) + Qctrl_ki = 3 / TimeConstant_Q; // Integral Gain of Q control + //+Vac control + X_tot = X_PCC_pu + Xtransfo_pu; // Total impedance in pu at wind generator connection point + VacCtrl_ki = 3 / (X_tot * TimeConstant_Vac); // Integral Gain of Vac control + VacCtrl_kp = 1; // Proportional Gain of Vac control (Only Integral part of the PI control is used) + //+P control + Pctrl_kp = 1; // Proportional Gain of P control (Only Integral part of the PI control is used) + Pctrl_ki = 3 / TimeConstant_P; // Integral Gain of P control + //+Maximum Power Point Tracking (MPPT) + Lambda_optimal = 6.635; // Optimal Lambda in order to have the rated active power for a wind speed = 11.24 m/s + Cp_max = 0.46; // Related Cp_max in order to have the rated active power for a wind speed = 11.24 m/s + K_optimal_pu = K_air * Cp_max * pow(K_speed / Lambda_optimal, 3) / Sgen_VA; // Constant coefficient defined as: P = K_optimal_pu * w^3 + //+ Idq limit (see CONTROL/Rotor_Control/Idq_ref_limiter1 ) + Q_priority = 1; // 1 = priority is given for reactive power (should be set for FRT capability). 0 = priority is given to active power + I_lim_Rotor_pu = 1.1; // Maximum rotor Idq limit + Id_lim_Rotor_pu = 1; // Maximum rotor Id limit, if Q_priority = 1 + Iq_lim_Rotor_pu = 1.1; // Maximum rotor Iq limit, if Q_priority = 0 + + //*Grid control (see CONTROL/Grid_Control ) + //+Grid Inner current control + eps = 0.7; // Damping ratio + w_i = 3 / (eps * TimeConstant_i_grid); // Impulse of the closed loop control + GridCtrl_ki = (Lchoke_pu / w) * w_i * w_i; // Integral Gain of the Inner Current control + GridCtrl_kp = 2 * eps * w_i * (Lchoke_pu / w); // Proportional gain of the Inner Current control + //+Vdc control + E_Cdc_J = 0.5 * CdcPark_F * Vdc_V * Vdc_V; // Energy stored in the DC capacitor (Joule) + H_Cdc = E_Cdc_J / (0.3 * Spark_VA); // Static moment of inertia (S). Note: Power converter is rated for 30% of the nominal power + eps = 1; // Damping ratio + w_Vdc = 3 / (eps * TimeConstant_Vdc); // Time response for 5% (sec) + VdcCtrl_ki = H_Cdc * w_Vdc * w_Vdc; // Proportional Gain of Vdc control + VdcCtrl_kp = 2 * eps * w_Vdc * H_Cdc; // Integral Gain of Vdc control + + //*PLL (see CONTROL/Computing_Variables/PLL) + PLL_lim_low = 12 + PLL_lim_up = 12 + PLL_max_delay = 0.0223 + eps = 1; // damping ratio + w_i = 3 / (eps * TimeConstant_PLL); + PLL_ki = 1 * w_i * w_i; // Int. Gain integral (rad/V) + PLL_kp = 2 * eps * w_i * 1; // Prop. Gain (rad/s/V) + + //*Protection system (see CONTROL/Protection ) + //+Deep Voltage Sag (DVS) protection + DVS_level = 0.001; // Threshold Voltage value in pu to activate the DVS protection + DVS_hysteresis = 0.1; // Threshold Voltage value in pu to deactivate the DVS protection + //+Crowbar protection + Rcrowbar_Ohm = 0.09; // Crowbar resistance (Ohm) + Lcrowbar_H = 4.0e-05; // Crowbar inductance (H) + RcrowbarPark_Ohm = Rcrowbar_Ohm / Ngen; // Crowbar equivalent resistance (Ohm) + LcrowbarPark_H = Lcrowbar_H / Ngen; // Crowbar equivalent inductance (H) + //+Converter Overcurrent limit + I_RotorSideConv_max = 4; // Threshold Current value in pu to activate RotorSide Overcurrent protection + I_GridSideConv_max = 4; // Threshold Current value in pu to activate GridSide Overcurrent protection + + //////////////////////////////////////////////////_INITIALIZATION ////////////////////////////////////////////////// + //*Initialization of the Protection system (see CONTROL/Protection) + init_Prot_delay = 0.3; // Delay time initialization of the Protection devices (s) (All protections, expect the Chopper, are deactivated for the specified time ) + + //*Initialization of the WIND TURBINE model (see WIND_TURBINE/Two_mass_model) + Lambda_ini_guess = Lambda_optimal; // Initial Lambda guess value + Cp_init_guess = Cp_max; // Initial Cp guess value + max_wind_speed = K_speed * w_rotor_max_pu / Lambda_ini_guess; // maximum wind speed before the Pitch control is activated + if (Mean_wind_speed > max_wind_speed) { // if the wind speed value specified by the user is higher than the maximum wind speed before Pitch control is activated + Mean_wind_speed_init = max_wind_speed; // than the wind speed value is fixed at his maximum value for initialization. + } + else { + Mean_wind_speed_init = Mean_wind_speed; + } + + init_w_turbine_pu = Mean_wind_speed_init * Lambda_ini_guess / K_speed * 0.989; // Initial Turbine speed in pu. The value is multiplied by a factor of losses. + init_Pturbine_pu = K_air * Cp_init_guess * pow(Mean_wind_speed_init, 3) / Sgen_VA; + init_Tturbine_pu = init_Pturbine_pu / init_w_turbine_pu; // Initial Turbine Torque in pu. + + //*Initialization of _ASYNCHRONOUS GENERATOR model (See ASM_DFIG/ASM1) + init_slip_pct = (1 - init_w_turbine_pu / 1) * 100; // Initial slip (%) + + /////////////////////////////////////////////////////////////_LOAD-FLOW CALCULATIONS ////////////////////////////////////////////////// + + //* Load-Flow of the wind generators (see LF_init1) + LF_V = Vgen_kVRMSLL * 1000 * sqrt(2 / 3); // Voltage for the LF device at the bus of the generator + + if (ControlMode == 1 || ControlMode == 3) { // if Q control or Q+Vac Control is selected, than PQ constraint is activated + LF_t1 = -1; // switch trick to active PQ constraint + LF_t2 = 1e15; // switch trick to active PQ constraint + LF_Ppark1_W = init_Pturbine_pu * Spark_VA * 0.967; // Active power in W of PQ constraint. The value is multiplied by a factor of losses. + LF_Ppark2_W = 0; // Active power in W of PV constraint. The value is set to zero + } else { // if Vac control is selected, than PV constraint is activated + LF_t1 = 1e15; // switch trick to active PV constraint + LF_t2 = -1; // switch trick to active PV constraint + LF_Ppark1_W = 0; // Active power in W of PQ constraint. The value is set to zero + LF_Ppark2_W = init_Pturbine_pu * Spark_VA * 0.967; // Active power in W of PV constraint. The value is multiplied by a factor of losses. + } + + +} + diff --git a/EMTP-RV/ADN EPEC2017/PV_maskver3.dwj b/EMTP-RV/ADN EPEC2017/PV_maskver3.dwj new file mode 100644 index 0000000..909e122 --- /dev/null +++ b/EMTP-RV/ADN EPEC2017/PV_maskver3.dwj @@ -0,0 +1,202 @@ +// Last change: +//*This script is used for the Photovoltaic (PV) model. +// This scripts allows to modify parameters that are not accessible on the main mask. + +with (Math) { + + //////////////////////////////////////////////////_BASIC CALCULATIONS ////////////////////////////////////////////////// + //* Fundamental Period + Period = 1 / Freq; // Period in s + period = Period; + w = 2 * PI * Freq; + //* Total solar panel S and P + Sgen_MVA = Pgen_MW * 1.11; // Rated apparent power of solar generator (MVA). Note that 11% of reserve power is added + Pgen_W = Pgen_MW * 1e6; // Rated active power of solar generator (W) + Sgen_VA = Sgen_MVA * 1e6; // Rated apparent power of solar generator (VA). Note that 11% of reserve power is added. + Spark_VA = Ngen * Sgen_VA; // Total aggregated solar plant S (VA) + Spark_MVA = Ngen * Sgen_VA / 1e6; // Total aggregated solar plant S (MVA), + Ppark_W = Ngen * Pgen_W; // Total aggregated solar plant P (W) (used for the load-flow device) + Qpark_Var = Qref_pu * Spark_VA; // Total aggregated solar plant Q (Var) (used for the load-flow device when Q or Q+Vac control is selected) + Vdc_V = Vdc_kV * 1e3; // DC voltage (V) + //* base calculations + Vbase_VRMSLL = Vgen_kVRMSLL * 1e3; // Rated generator voltage (VRMSLL) + Vbase_V = Vbase_VRMSLL * sqrt(2 / 3); // Peak generator voltage (V) + Ibase_A = Spark_VA / Vbase_VRMSLL * sqrt(2 / 3); // Peak generator current (A) + Zbase_Ohm = (Vgen_kVRMSLL * Vgen_kVRMSLL) / Spark_MVA; // Base impedance (Ohm) + + //////////////////////////////////////////////////_ DC/AC CONVERTER PARAMETERS ////////////////////////////////////////////////// + //* Choke impedance (see RLchoke device) + Rchoke_pu = 1.5e-3; // in pu + Rchoke_Ohm = Rchoke_pu * Zbase_Ohm; // in Ohm + Lchoke_pu = 0.15; // in pu + Lchoke_H = Lchoke_pu * Zbase_Ohm / w; // in H + + //* DC/AC Converter parameters (see Converter device) + E_Cdc = 10; // in KJ/MVA + Cdc_F = E_Cdc * 1e-3 / (1 / 2 * pow(Vdc_V, 2) / Sgen_VA); // DC capacitor (F) + CdcPark_F = Cdc_F * Ngen; // Aggregated DC capacitor (F) + CarrierSignal_Freq = 3e3; // PWM Carrier signal frequency (Hz) (used only for Detailed Converter Model) + CarrierSignal_ratio = CarrierSignal_Freq / Freq; // PWM carrier signal ratio (used only for Detailed Converter Model) + Rsnubber_Ohm = 3000; // IGBT/diode resistance snubber (Ohm) (used only for Detailed Converter Model) + Csnubber_F = 1e-6; // IGBT/diode capacitor snubber (F) (used only for Detailed Converter Model) + Rchopper_Ohm = 0.2; // Chopper resistance in Ohm + RchopperPark_Ohm = Rchopper_Ohm / Ngen; // Aggregation of the chopper resistance + + //* High pass filter parameters (see filter1 device) + Qfilter_Var = 150e3 * Ngen // Reactive power generated by each filter + n1 = CarrierSignal_ratio // First tuning harmonic + n2 = CarrierSignal_ratio * 2 // Second tuning harmonic + QF = 5 // Quality factor + Cfilter1_F = Qfilter_Var / (Vbase_VRMSLL * Vbase_VRMSLL * 2 * PI * Freq); // capacitor in F + Lfilter1_H = 1 / (Cfilter1_F * pow((n1 * Freq * 2 * PI), 2)); // inductance in H + Rfilter1_Ohm = QF * 2 * PI * n1 * Freq * Lfilter1_H; // resistance in Ohm + Cfilter2_F = Cfilter1_F; // capacitor in F + Lfilter2_H = 1 / (Cfilter2_F * pow((n2 * Freq * 2 * PI), 2)); // inductance in H + Rfilter2_Ohm = QF * 2 * PI * n2 * Freq * Lfilter2_H; // resistance in Ohm + + + //////////////////////////////////////////////////SOLAR INVERTER AND _TRANSFORMER PARAMETERS ////////////////////////////////////////////////// + //* solar inverter transformer parameters (see Transformer device) + Rtransfo_sec_Ohm = (Rtransfo_pu * 0.1 * Vgrid_kVRMSLL * Vgrid_kVRMSLL / Spark_MVA); // Secondary winding resistance, Ohm + Ltransfo_sec_H = (Xtransfo_pu * 0.1 * Vgrid_kVRMSLL * Vgrid_kVRMSLL / Spark_MVA) / (2 * PI * Freq); // Secondary winding inductance, H + Rtransfo_prim_Ohm = (Rtransfo_pu * 0.9 * Vgen_kVRMSLL * Vgen_kVRMSLL / Spark_MVA); // Primary winding resistance, Ohm + Ltransfo_prim_H = (Xtransfo_pu * 0.9 * Vgen_kVRMSLL * Vgen_kVRMSLL / Spark_MVA) / (2 * PI * Freq); // Primary winding inductance, H + Transfo_ratio = Vgen_kVRMSLL / Vgrid_kVRMSLL / sqrt(3); // Turns ratio + + + //////////////////////////////////////////////////_CONTROL SYSTEM PARAMETERS ////////////////////////////////////////////////// + //* Grid control (see CONTROL/Grid_Control ) + //+Grid Inner current control + eps = 1; // Damping ratio + w_i = 3 / (eps * TimeConstant_i_grid); // Impulse of the closed loop control + GridCtrl_ki = (Lchoke_pu / w) * w_i * w_i; // Integral Gain of the Inner Current control + GridCtrl_kp = 5 * 2 * eps * w_i * (Lchoke_pu / w); // Proportional gain of the Inner Current control + //+Vdc control + E_Cdc_J = 0.5 * CdcPark_F * Vdc_V * Vdc_V; // Energy stored in the DC capacitor (Joule) + H_Cdc = E_Cdc_J / (1 * Ppark_W); // Static moment of inertia (S) + eps = 0.7; // Damping ratio + w_Vdc = 3 / (eps * TimeConstant_Vdc); // Time response for 5% (sec) + VdcCtrl_ki = 2 * H_Cdc * w_Vdc * w_Vdc; // Proportional Gain of Vdc control + VdcCtrl_kp = 2 * 2 * eps * w_Vdc * H_Cdc; // Integral Gain of Vdc control + //+ Idq limit (see CONTROL/Grid_Control/Idq_ref_limiter1 ) + Q_priority = 1; // 1 = priority is given for reactive power (should be set for FRT capability). 0 = priority is given to active power + I_lim_Grid_pu = 1.1; // Maximum Grid Converter Idq limit + Id_lim_Grid_pu = 1; // Maximum Grid Converter Id limit, if Q_priority = 1 + Iq_lim_Grid_pu = 1; // Maximum Grid Converter Iq limit, if Q_priority = 0 + //+Q control and Q+Vac control + Qctrl_kp = 0; // Proportional Gain of Q control (Only Integral part of the PI control is used) + Qctrl_ki = 3 / TimeConstant_Q; // Integral Gain of Q control + //+Vac control + X_tot = X_PCC_pu + Xtransfo_pu; // Total impedance in pu at solar generator connection point + VacCtrl_kp = 1; // Proportional Gain of Vac control (Only Integral part of the PI control is used) + VacCtrl_ki = 3 / (X_tot * TimeConstant_Vac); // Integral Gain of Vac control + + //* PLL (see CONTROL/Computing_Variables/PLL) + PLL_lim_low = 12 + PLL_lim_up = 12 + PLL_max_delay = 0.0223 + eps = 1; // Damping ratio + w_i = 3 / (eps * TimeConstant_PLL); + PLL_ki = 1 * w_i * w_i; // Integral Gain (rad/V) + PLL_kp = 2 * eps * w_i * 1; // Proportional Gain (rad/s/V) + + //* Protection system (see CONTROL/Protection ) + //+Deep Voltage Sag (DVS) protection + DVS_level = 0.001; // Threshold Voltage value in pu to activate the DVS protection + DVS_hysteresis = 0.1; // Threshold Voltage value in pu to deactivate the DVS protection + //+Converter Overcurrent limit + I_GridSideConv_max = 4; // Threshold Current value in pu to activate GridSide Overcurrent protection + + //////////////////////////////////////////////////_INITIALIZATION ////////////////////////////////////////////////// + //* Initialization of the Protection system (see CONTROL/Protection) + init_Prot_delay = 0.3; // Delay time initialization of the Protection devices (s) (All protections, expect the Chopper, are deactivated for the specified time ) + + /////////////////////////////////////////////////////////////_LOAD-FLOW CALCULATIONS ////////////////////////////////////////////////// + //* Load-Flow of the solar plant (see LF_init1) + LF_V = Vref_pu * Vgen_kVRMSLL * 1000 * sqrt(2 / 3); // Voltage for the LF device at the bus of the generator + + if (ControlMode == 1 || ControlMode == 3) { // if Q control or Q+Vac Control is selected, than PQ constraint is activated + LF_t1 = -1; // switch trick to active PQ constraint + LF_t2 = 1e15; // switch trick to active PQ constraint + LF_Ppark1_W = Spark_VA / (1.1162); // Active power in W of PQ constraint. The value is multiplied to account for losses and to remove 1.11 reserve factor also + LF_Ppark2_W = 0; // Active power in W of PV constraint. The value is set to zero + } else { // if Vac control is selected, than PV constraint is activated + LF_t1 = 1e15; // switch trick to active PV constraint + LF_t2 = -1; // switch trick to active PV constraint + LF_Ppark1_W = 0; // Active power in W of PQ constraint. The value is set to zero + LF_Ppark2_W = Spark_VA / (1.1162); // Active power in W of PV constraint. The value is multiplied to account for losses and to remove 1.11 reserve factor also + } + + + + + + // Q/V Regulator _ modified controller (updated version of VQ controller to calculate Iq_ref) + kp_var = 0.5; //0 + ki_var = 60; //40// + kp_volt = 1 // 0 + ki_volt = 60 // + + VoltReg_FRT_Gain = 2; // Proportional Gain of V control in FRT function + FRT_ON = 0.1; // FRT function is activated when voltage deviation from 1 pu (dV) is larger than FRT_ON + FRT_OFF = 0.075; // FRT function is deactivated when (dV < FRT_OFF) & (t_FRT > FRT_time) + FRT_time = 0.25; + + + + + + ////////////////////////////////////////////////////SOLAR PANEL PARAMETERS ////////////////////////////////// + /// ~75 W per module + // input + // Calcul + + //dVmpp=16.2; // maximum power voltage (V) under STC + //dnNbPVCellsSeries=Vdc_V / dVmpp; + //PoweratVmpp=73.9273; // new model diode + //dnNbPVCellsParallel=Pgen_W/dnNbPVCellsSeries/PoweratVmpp; + + Nmod_series = Vdc_V / VmaxP; // number of PV modules in series + Nmod_parallel = Ppark_W/(Vdc_V*ImaxP); // number of PV modules in parallel + //Nmod_parallel = Spark_VA/(Vdc_V*ImaxP); // number of PV modules in parallel + + ///////////// SOLAR PANEL PARAMETER CALCULATION from DATASHEET /////// (13/02/2015 Thomas Kauffmann) + + Vth = (Temp_ref + 273)*1.38065e-23/1.6022e-19; // diode threshold voltage + Ns = Ncell_series; + Io = Isc / (-1 + exp(Voc /(IdealFactor*Ns*Vth))); // diode reversed saturation current + Rs_max = Ns*((Voc - VmaxP)/(Ns*ImaxP) - exp(-Voc/(Vth*Ns*IdealFactor))*Vth*IdealFactor/Io); + + x1 = 10*Rs_max; // to initialize the while loop + x2 = Rs_max; + ite = 0; + + while ((x1 - x2)>num_tolerance){ + x1 = x2; + Rp_x1 = 1/(1/(VmaxP/ImaxP - x1) - Io/(IdealFactor*Ns*Vth)*exp((VmaxP + ImaxP*x1)/(IdealFactor*Ns*Vth))); + Iph_x1 = Isc*(VmaxP/(VmaxP - x1*ImaxP) - Io*x1*exp((VmaxP + ImaxP*x1)/(IdealFactor*Ns*Vth))/IdealFactor/Ns/Vth); + + // Newton Method + fx1 = Iph_x1 - ImaxP - (VmaxP + ImaxP*x1)/Rp_x1 - Io*(exp((VmaxP + ImaxP*x1)/(IdealFactor*Ns*Vth))-1); + df_x1 = ImaxP*VmaxP*(Isc - 2*ImaxP)/pow(VmaxP - x1*ImaxP,2) + Io*exp((VmaxP + ImaxP*x1)/(IdealFactor*Ns*Vth))*(ImaxP*VmaxP - Isc*IdealFactor*Ns*Vth + x1*ImaxP*(ImaxP - Isc))/pow(IdealFactor*Ns*Vth,2); + x2 = x1 - fx1/df_x1; + + ite = ite+1; + + } + + Rs = x2; + Rp = 1/(1/(VmaxP/ImaxP - x2) - Io/(IdealFactor*Ns*Vth)*exp((VmaxP + ImaxP*x2)/(IdealFactor*Ns*Vth))); + Iph = Isc*(Rs + Rp)/Rp; + + Rseries_PV = Rs*(Nmod_series/Nmod_parallel); + Rparallel_PV = Rp*(Nmod_series/Nmod_parallel); + + Vth_diode = (Temp + 273)*1.38065e-23/1.6022e-19; + I0_diode = Nmod_parallel*(Isc + Ki_pv*(Temp - Temp_ref))/(-1 + exp((Voc + Kv_pv*(Temp - Temp_ref))/(IdealFactor*Ns*Vth))); + Iph_T = Nmod_parallel*(Iph + Ki_pv*(Temp - Temp_ref)); + + Nseries = Ns*Nmod_series // total number of cell in series (value transmitted) + + +} diff --git a/EMTP-RV/ADN EPEC2017/kth_EPEC2017.ecf b/EMTP-RV/ADN EPEC2017/kth_EPEC2017.ecf new file mode 100644 index 0000000..ebc6093 Binary files /dev/null and b/EMTP-RV/ADN EPEC2017/kth_EPEC2017.ecf differ diff --git a/EMTP-RV/ADN IPST2017/DFIG_mask.dwj b/EMTP-RV/ADN IPST2017/DFIG_mask.dwj new file mode 100644 index 0000000..7709f84 --- /dev/null +++ b/EMTP-RV/ADN IPST2017/DFIG_mask.dwj @@ -0,0 +1,201 @@ +// Last change: HANI SAAD - JEAN MAHSEREDJIAN 04/11/2013 03:36:19 AM +//*This script is used for the DFIG wind generator model. +// This scripts allows to modify parameters that are not accessible on the main mask. + +with (Math) { + + //////////////////////////////////////////////////_BASIC CALCULATIONS ////////////////////////////////////////////////// + //*Fundamental Period + Period = 1 / Freq; + period = Period; + w = 2 * PI * Freq; + //*Total Wind park S and P + Pgen_W = Pgen_MW * 1e6; // Rated active power of wind generator (W) + Sgen_VA = Pgen_W * 1.11; // Rated apparent power of wind generator (VA). Note that 11% of reserve power is added. + Spark_VA = Ngen * Sgen_VA; // Total aggregated wind park S (VA) + Spark_MVA = Ngen * Sgen_VA / 1e6; // Total aggregated wind park S (MVA), ASM machine uses MVA selection + Ppark_W = Ngen * Pgen_W; // Total aggregated wind park P (W) (used for the load-flow device) + Qpark_Var = Qref_pu * Spark_VA; // Total aggregated wind park Q (Var) (used for the load-flow device when Q or Q+Vac control is selected) + Vdc_V = Vdc_kV * 1e3; // DC voltage (V) + //*base calculations + Vbase_VRMSLL = Vgen_kVRMSLL * 1e3; // Rated generator voltage (VRMSLL) + Vbase_V = Vbase_VRMSLL * sqrt(2 / 3); // Peak generator voltage (V) + Ibase_A = Spark_VA / Vbase_VRMSLL * sqrt(2 / 3); // Peak generator current (A) + Zbase_Ohm = (Vgen_kVRMSLL * Vgen_kVRMSLL) / Spark_MVA; // Base impedance (Ohm) + + //////////////////////////////////////////////////_ASYNCHRONOUS GENERATOR & AC/DC/AC CONVERTER PARAMETERS ////////////////////////////////////////////////// + //*ASM parameters (see ASM_DFIG/ASM1) + npoles = 6; // Number of poles of the asynchronous generator + Rs_pu = 0.007 // in pu + Lls_pu = 0.171; // in pu + Lmd_pu = 2.9 // in pu + Lmq_pu = 2.9 // in pu + Rr_pu = 0.005 // in pu + Llr_pu = 0.156 // in pu + Hgen_pu = 5.04 // in sec + D_pu = 0.0 // in pu + Llr_plus_Lmd_pu = Llr_pu + Lmd_pu; // in pu (used in DFIG_CONTROL/ROTOR side CONTROL) + //*Choke impedance (see ASM_DFIG/Choke) + Rchoke_pu = 0.003; // in pu + Rchoke_Ohm = Rchoke_pu * Zbase_Ohm; // in Ohm + Lchoke_pu = 0.3 // in pu + Lchoke_H = Lchoke_pu * Zbase_Ohm / (2 * PI * Freq); // in H + //*AC/DC/AC Converter parameters (see ASM_DFIG/BtoB_converter) + E_Cdc = 5; // in KJ/MVA + Cdc_F = E_Cdc * 1e-3 / (1 / 2 * pow(Vdc_V, 2) / Sgen_VA); // DC capacitor (F) + CdcPark_F = Cdc_F * Ngen; // Aggregated DC capacitor (F) + CarrierSignal_Freq = 3e3; // PWM Carrier signal frequency (Hz) (used only for Detailed Converter Model) + CarrierSignal_ratio = CarrierSignal_Freq / Freq; // PWM carrier signal ratio (used only for Detailed Converter Model) + Rsnubber_Ohm = 3000; // IGBT/diode resistance snubber (Ohm) (used only for Detailed Converter Model) + Csnubber_F = 1e-6; // IGBT/diode capacitor snubber (F) (used only for Detailed Converter Model) + Rchopper_Ohm = 0.2; // Chopper resistance in Ohm + RchopperPark_Ohm = Rchopper_Ohm / Ngen; // Aggregation of the chopper resistance + //*High pass filter parameters (see ASM_DFIG/filter_stator_ASM ) + Qfilter_Var = 150e3 * Ngen // Reactive power generated by each filter + n1 = CarrierSignal_ratio // First tuning harmonic + n2 = CarrierSignal_ratio * 2 // Second tuning harmonic + QF = 1e3 // Quality factor + Cfilter1_F = Qfilter_Var / (Vbase_VRMSLL * Vbase_VRMSLL * 2 * PI * Freq); + Lfilter1_H = 1 / (Cfilter1_F * pow((n1 * Freq * 2 * PI), 2)); + Rfilter1_Ohm = QF * 2 * PI * n1 * Freq * Lfilter1_H; + Cfilter2_F = Cfilter1_F; + Lfilter2_H = 1 / (Cfilter2_F * pow((n2 * Freq * 2 * PI), 2)); + Rfilter2_Ohm = QF * 2 * PI * n2 * Freq * Lfilter2_H; + + //////////////////////////////////////////////////_WIND TURBINE PARAMETERS ////////////////////////////////////////////////// + //* Blade and Gearbox parameters (see WIND_TURBINE/Wind_Power) + blade_radius = ceil((sqrt(pow(35.25, 2) / 1.5e6 * Pgen_W)) * 100) / 100; // Blade radius (m) is in function of the rated power + gearbox_ratio = 72 * (Freq / 60); // Gearbox ratio. Note that the ratio is adjusted depending on the frequency used in order to preserve the same tip-speed value (Lambda). + ro = 1.225; // air density (kg/m^3) + K_speed = 2 * PI * Freq / (npoles / 2) * blade_radius / gearbox_ratio; // speed coefficient defined as: Lambda = K_speed / wind_speed * w_rotor; + K_air = (1 / 2) * PI * blade_radius * blade_radius * ro; // air coefficient defined as: Pair = K_air * Cp * (wind_speed)^3; + Pitch_max_deg = 20; // Maximum pitch angle in degree + Pitch_min_deg = -1.5; // Minimum pitch angle in degree + + //*Shaft and train drive parameters (see WIND_TURBINE/Two_mass_model) + Ksh = 1.2; //Speed deviation damping + Dm = 1.5; //Mutual damping + Ht = 4.0 // Inertia constant (s) + + //////////////////////////////////////////////////_TRANSFORMER PARAMETERS ////////////////////////////////////////////////// + //* Wind generator transformer parameters (see Transformer) + Rtransfo_sec_Ohm = (Rtransfo_pu * 0.1 * Vgrid_kVRMSLL * Vgrid_kVRMSLL / Spark_MVA); // Secondary winding resistance, Ohm + Ltransfo_sec_H = (Xtransfo_pu * 0.1 * Vgrid_kVRMSLL * Vgrid_kVRMSLL / Spark_MVA) / (2 * PI * Freq); // Secondary winding inductance, H + Rtransfo_prim_Ohm = (Rtransfo_pu * 0.9 * Vgen_kVRMSLL * Vgen_kVRMSLL / Spark_MVA); // Primary winding resistance, Ohm + Ltransfo_prim_H = (Xtransfo_pu * 0.9 * Vgen_kVRMSLL * Vgen_kVRMSLL / Spark_MVA) / (2 * PI * Freq); // Primary winding inductance, H + Transfo_ratio = Vgen_kVRMSLL / Vgrid_kVRMSLL / sqrt(3); // Turns ratio + + //////////////////////////////////////////////////_CONTROL SYSTEM PARAMETERS ////////////////////////////////////////////////// + //* Pitch control (see CONTROL/Pitch_Control ) + // NB: Ptich Control Gains are tuned by trial and error + PitchCtrl_ki = 20; // Integral Gain of the Pitch control + PitchCtrl_kp = 50; // Proportional Gain of the Pitch control + PitchCtrl_Td_d = 10; // Derivative Gain of the Pitch control + PitchCtrl_tau_d = 5e-3; // derivative filter in s + w_rotor_max_pu = 1.25; // maximum rotor speed in pu + Pitch_rate = 10 // Slop rate (see Pitch_Control/Rate_Limiter) + + //*Rotor control (see CONTROL/Rotor_Control ) + //+Rotor Inner current control + eps = 0.7; // Damping ratio + w_i = 3 / (eps * TimeConstant_i_rotor); // Impulse of the closed loop control + RotorCtrl_ki = 0.2 * ((Llr_plus_Lmd_pu) / w) * w_i * w_i; // Integral Gain of the Inner Current control + RotorCtrl_kp = 0.2 * 2 * eps * w_i * ((Llr_plus_Lmd_pu) / w); // Proportional gain of the Inner Current control + //+Q control and Q+Vac control + Qctrl_kp = 0; // Proportional Gain of Q control (Only Integral part of the PI control is used) + Qctrl_ki = 3 / TimeConstant_Q; // Integral Gain of Q control + //+Vac control + X_tot = X_PCC_pu + Xtransfo_pu; // Total impedance in pu at wind generator connection point + VacCtrl_ki = 3 / (X_tot * TimeConstant_Vac); // Integral Gain of Vac control + VacCtrl_kp = 1; // Proportional Gain of Vac control (Only Integral part of the PI control is used) + //+P control + Pctrl_kp = 1; // Proportional Gain of P control (Only Integral part of the PI control is used) + Pctrl_ki = 3 / TimeConstant_P; // Integral Gain of P control + //+Maximum Power Point Tracking (MPPT) + Lambda_optimal = 6.635; // Optimal Lambda in order to have the rated active power for a wind speed = 11.24 m/s + Cp_max = 0.46; // Related Cp_max in order to have the rated active power for a wind speed = 11.24 m/s + K_optimal_pu = K_air * Cp_max * pow(K_speed / Lambda_optimal, 3) / Sgen_VA; // Constant coefficient defined as: P = K_optimal_pu * w^3 + //+ Idq limit (see CONTROL/Rotor_Control/Idq_ref_limiter1 ) + Q_priority = 1; // 1 = priority is given for reactive power (should be set for FRT capability). 0 = priority is given to active power + I_lim_Rotor_pu = 1.1; // Maximum rotor Idq limit + Id_lim_Rotor_pu = 1; // Maximum rotor Id limit, if Q_priority = 1 + Iq_lim_Rotor_pu = 1.1; // Maximum rotor Iq limit, if Q_priority = 0 + + //*Grid control (see CONTROL/Grid_Control ) + //+Grid Inner current control + eps = 0.7; // Damping ratio + w_i = 3 / (eps * TimeConstant_i_grid); // Impulse of the closed loop control + GridCtrl_ki = (Lchoke_pu / w) * w_i * w_i; // Integral Gain of the Inner Current control + GridCtrl_kp = 2 * eps * w_i * (Lchoke_pu / w); // Proportional gain of the Inner Current control + //+Vdc control + E_Cdc_J = 0.5 * CdcPark_F * Vdc_V * Vdc_V; // Energy stored in the DC capacitor (Joule) + H_Cdc = E_Cdc_J / (0.3 * Spark_VA); // Static moment of inertia (S). Note: Power converter is rated for 30% of the nominal power + eps = 1; // Damping ratio + w_Vdc = 3 / (eps * TimeConstant_Vdc); // Time response for 5% (sec) + VdcCtrl_ki = H_Cdc * w_Vdc * w_Vdc; // Proportional Gain of Vdc control + VdcCtrl_kp = 2 * eps * w_Vdc * H_Cdc; // Integral Gain of Vdc control + + //*PLL (see CONTROL/Computing_Variables/PLL) + PLL_lim_low = 12 + PLL_lim_up = 12 + PLL_max_delay = 0.0223 + eps = 1; // damping ratio + w_i = 3 / (eps * TimeConstant_PLL); + PLL_ki = 1 * w_i * w_i; // Int. Gain integral (rad/V) + PLL_kp = 2 * eps * w_i * 1; // Prop. Gain (rad/s/V) + + //*Protection system (see CONTROL/Protection ) + //+Deep Voltage Sag (DVS) protection + DVS_level = 0.001; // Threshold Voltage value in pu to activate the DVS protection + DVS_hysteresis = 0.1; // Threshold Voltage value in pu to deactivate the DVS protection + //+Crowbar protection + Rcrowbar_Ohm = 0.09; // Crowbar resistance (Ohm) + Lcrowbar_H = 4.0e-05; // Crowbar inductance (H) + RcrowbarPark_Ohm = Rcrowbar_Ohm / Ngen; // Crowbar equivalent resistance (Ohm) + LcrowbarPark_H = Lcrowbar_H / Ngen; // Crowbar equivalent inductance (H) + //+Converter Overcurrent limit + I_RotorSideConv_max = 4; // Threshold Current value in pu to activate RotorSide Overcurrent protection + I_GridSideConv_max = 4; // Threshold Current value in pu to activate GridSide Overcurrent protection + + //////////////////////////////////////////////////_INITIALIZATION ////////////////////////////////////////////////// + //*Initialization of the Protection system (see CONTROL/Protection) + init_Prot_delay = 0.3; // Delay time initialization of the Protection devices (s) (All protections, expect the Chopper, are deactivated for the specified time ) + + //*Initialization of the WIND TURBINE model (see WIND_TURBINE/Two_mass_model) + Lambda_ini_guess = Lambda_optimal; // Initial Lambda guess value + Cp_init_guess = Cp_max; // Initial Cp guess value + max_wind_speed = K_speed * w_rotor_max_pu / Lambda_ini_guess; // maximum wind speed before the Pitch control is activated + if (Mean_wind_speed > max_wind_speed) { // if the wind speed value specified by the user is higher than the maximum wind speed before Pitch control is activated + Mean_wind_speed_init = max_wind_speed; // than the wind speed value is fixed at his maximum value for initialization. + } + else { + Mean_wind_speed_init = Mean_wind_speed; + } + + init_w_turbine_pu = Mean_wind_speed_init * Lambda_ini_guess / K_speed * 0.989; // Initial Turbine speed in pu. The value is multiplied by a factor of losses. + init_Pturbine_pu = K_air * Cp_init_guess * pow(Mean_wind_speed_init, 3) / Sgen_VA; + init_Tturbine_pu = init_Pturbine_pu / init_w_turbine_pu; // Initial Turbine Torque in pu. + + //*Initialization of _ASYNCHRONOUS GENERATOR model (See ASM_DFIG/ASM1) + init_slip_pct = (1 - init_w_turbine_pu / 1) * 100; // Initial slip (%) + + /////////////////////////////////////////////////////////////_LOAD-FLOW CALCULATIONS ////////////////////////////////////////////////// + + //* Load-Flow of the wind generators (see LF_init1) + LF_V = Vgen_kVRMSLL * 1000 * sqrt(2 / 3); // Voltage for the LF device at the bus of the generator + + if (ControlMode == 1 || ControlMode == 3) { // if Q control or Q+Vac Control is selected, than PQ constraint is activated + LF_t1 = -1; // switch trick to active PQ constraint + LF_t2 = 1e15; // switch trick to active PQ constraint + LF_Ppark1_W = init_Pturbine_pu * Spark_VA * 0.967; // Active power in W of PQ constraint. The value is multiplied by a factor of losses. + LF_Ppark2_W = 0; // Active power in W of PV constraint. The value is set to zero + } else { // if Vac control is selected, than PV constraint is activated + LF_t1 = 1e15; // switch trick to active PV constraint + LF_t2 = -1; // switch trick to active PV constraint + LF_Ppark1_W = 0; // Active power in W of PQ constraint. The value is set to zero + LF_Ppark2_W = init_Pturbine_pu * Spark_VA * 0.967; // Active power in W of PV constraint. The value is multiplied by a factor of losses. + } + + +} + diff --git a/EMTP-RV/ADN IPST2017/PV_maskver3.dwj b/EMTP-RV/ADN IPST2017/PV_maskver3.dwj new file mode 100644 index 0000000..909e122 --- /dev/null +++ b/EMTP-RV/ADN IPST2017/PV_maskver3.dwj @@ -0,0 +1,202 @@ +// Last change: +//*This script is used for the Photovoltaic (PV) model. +// This scripts allows to modify parameters that are not accessible on the main mask. + +with (Math) { + + //////////////////////////////////////////////////_BASIC CALCULATIONS ////////////////////////////////////////////////// + //* Fundamental Period + Period = 1 / Freq; // Period in s + period = Period; + w = 2 * PI * Freq; + //* Total solar panel S and P + Sgen_MVA = Pgen_MW * 1.11; // Rated apparent power of solar generator (MVA). Note that 11% of reserve power is added + Pgen_W = Pgen_MW * 1e6; // Rated active power of solar generator (W) + Sgen_VA = Sgen_MVA * 1e6; // Rated apparent power of solar generator (VA). Note that 11% of reserve power is added. + Spark_VA = Ngen * Sgen_VA; // Total aggregated solar plant S (VA) + Spark_MVA = Ngen * Sgen_VA / 1e6; // Total aggregated solar plant S (MVA), + Ppark_W = Ngen * Pgen_W; // Total aggregated solar plant P (W) (used for the load-flow device) + Qpark_Var = Qref_pu * Spark_VA; // Total aggregated solar plant Q (Var) (used for the load-flow device when Q or Q+Vac control is selected) + Vdc_V = Vdc_kV * 1e3; // DC voltage (V) + //* base calculations + Vbase_VRMSLL = Vgen_kVRMSLL * 1e3; // Rated generator voltage (VRMSLL) + Vbase_V = Vbase_VRMSLL * sqrt(2 / 3); // Peak generator voltage (V) + Ibase_A = Spark_VA / Vbase_VRMSLL * sqrt(2 / 3); // Peak generator current (A) + Zbase_Ohm = (Vgen_kVRMSLL * Vgen_kVRMSLL) / Spark_MVA; // Base impedance (Ohm) + + //////////////////////////////////////////////////_ DC/AC CONVERTER PARAMETERS ////////////////////////////////////////////////// + //* Choke impedance (see RLchoke device) + Rchoke_pu = 1.5e-3; // in pu + Rchoke_Ohm = Rchoke_pu * Zbase_Ohm; // in Ohm + Lchoke_pu = 0.15; // in pu + Lchoke_H = Lchoke_pu * Zbase_Ohm / w; // in H + + //* DC/AC Converter parameters (see Converter device) + E_Cdc = 10; // in KJ/MVA + Cdc_F = E_Cdc * 1e-3 / (1 / 2 * pow(Vdc_V, 2) / Sgen_VA); // DC capacitor (F) + CdcPark_F = Cdc_F * Ngen; // Aggregated DC capacitor (F) + CarrierSignal_Freq = 3e3; // PWM Carrier signal frequency (Hz) (used only for Detailed Converter Model) + CarrierSignal_ratio = CarrierSignal_Freq / Freq; // PWM carrier signal ratio (used only for Detailed Converter Model) + Rsnubber_Ohm = 3000; // IGBT/diode resistance snubber (Ohm) (used only for Detailed Converter Model) + Csnubber_F = 1e-6; // IGBT/diode capacitor snubber (F) (used only for Detailed Converter Model) + Rchopper_Ohm = 0.2; // Chopper resistance in Ohm + RchopperPark_Ohm = Rchopper_Ohm / Ngen; // Aggregation of the chopper resistance + + //* High pass filter parameters (see filter1 device) + Qfilter_Var = 150e3 * Ngen // Reactive power generated by each filter + n1 = CarrierSignal_ratio // First tuning harmonic + n2 = CarrierSignal_ratio * 2 // Second tuning harmonic + QF = 5 // Quality factor + Cfilter1_F = Qfilter_Var / (Vbase_VRMSLL * Vbase_VRMSLL * 2 * PI * Freq); // capacitor in F + Lfilter1_H = 1 / (Cfilter1_F * pow((n1 * Freq * 2 * PI), 2)); // inductance in H + Rfilter1_Ohm = QF * 2 * PI * n1 * Freq * Lfilter1_H; // resistance in Ohm + Cfilter2_F = Cfilter1_F; // capacitor in F + Lfilter2_H = 1 / (Cfilter2_F * pow((n2 * Freq * 2 * PI), 2)); // inductance in H + Rfilter2_Ohm = QF * 2 * PI * n2 * Freq * Lfilter2_H; // resistance in Ohm + + + //////////////////////////////////////////////////SOLAR INVERTER AND _TRANSFORMER PARAMETERS ////////////////////////////////////////////////// + //* solar inverter transformer parameters (see Transformer device) + Rtransfo_sec_Ohm = (Rtransfo_pu * 0.1 * Vgrid_kVRMSLL * Vgrid_kVRMSLL / Spark_MVA); // Secondary winding resistance, Ohm + Ltransfo_sec_H = (Xtransfo_pu * 0.1 * Vgrid_kVRMSLL * Vgrid_kVRMSLL / Spark_MVA) / (2 * PI * Freq); // Secondary winding inductance, H + Rtransfo_prim_Ohm = (Rtransfo_pu * 0.9 * Vgen_kVRMSLL * Vgen_kVRMSLL / Spark_MVA); // Primary winding resistance, Ohm + Ltransfo_prim_H = (Xtransfo_pu * 0.9 * Vgen_kVRMSLL * Vgen_kVRMSLL / Spark_MVA) / (2 * PI * Freq); // Primary winding inductance, H + Transfo_ratio = Vgen_kVRMSLL / Vgrid_kVRMSLL / sqrt(3); // Turns ratio + + + //////////////////////////////////////////////////_CONTROL SYSTEM PARAMETERS ////////////////////////////////////////////////// + //* Grid control (see CONTROL/Grid_Control ) + //+Grid Inner current control + eps = 1; // Damping ratio + w_i = 3 / (eps * TimeConstant_i_grid); // Impulse of the closed loop control + GridCtrl_ki = (Lchoke_pu / w) * w_i * w_i; // Integral Gain of the Inner Current control + GridCtrl_kp = 5 * 2 * eps * w_i * (Lchoke_pu / w); // Proportional gain of the Inner Current control + //+Vdc control + E_Cdc_J = 0.5 * CdcPark_F * Vdc_V * Vdc_V; // Energy stored in the DC capacitor (Joule) + H_Cdc = E_Cdc_J / (1 * Ppark_W); // Static moment of inertia (S) + eps = 0.7; // Damping ratio + w_Vdc = 3 / (eps * TimeConstant_Vdc); // Time response for 5% (sec) + VdcCtrl_ki = 2 * H_Cdc * w_Vdc * w_Vdc; // Proportional Gain of Vdc control + VdcCtrl_kp = 2 * 2 * eps * w_Vdc * H_Cdc; // Integral Gain of Vdc control + //+ Idq limit (see CONTROL/Grid_Control/Idq_ref_limiter1 ) + Q_priority = 1; // 1 = priority is given for reactive power (should be set for FRT capability). 0 = priority is given to active power + I_lim_Grid_pu = 1.1; // Maximum Grid Converter Idq limit + Id_lim_Grid_pu = 1; // Maximum Grid Converter Id limit, if Q_priority = 1 + Iq_lim_Grid_pu = 1; // Maximum Grid Converter Iq limit, if Q_priority = 0 + //+Q control and Q+Vac control + Qctrl_kp = 0; // Proportional Gain of Q control (Only Integral part of the PI control is used) + Qctrl_ki = 3 / TimeConstant_Q; // Integral Gain of Q control + //+Vac control + X_tot = X_PCC_pu + Xtransfo_pu; // Total impedance in pu at solar generator connection point + VacCtrl_kp = 1; // Proportional Gain of Vac control (Only Integral part of the PI control is used) + VacCtrl_ki = 3 / (X_tot * TimeConstant_Vac); // Integral Gain of Vac control + + //* PLL (see CONTROL/Computing_Variables/PLL) + PLL_lim_low = 12 + PLL_lim_up = 12 + PLL_max_delay = 0.0223 + eps = 1; // Damping ratio + w_i = 3 / (eps * TimeConstant_PLL); + PLL_ki = 1 * w_i * w_i; // Integral Gain (rad/V) + PLL_kp = 2 * eps * w_i * 1; // Proportional Gain (rad/s/V) + + //* Protection system (see CONTROL/Protection ) + //+Deep Voltage Sag (DVS) protection + DVS_level = 0.001; // Threshold Voltage value in pu to activate the DVS protection + DVS_hysteresis = 0.1; // Threshold Voltage value in pu to deactivate the DVS protection + //+Converter Overcurrent limit + I_GridSideConv_max = 4; // Threshold Current value in pu to activate GridSide Overcurrent protection + + //////////////////////////////////////////////////_INITIALIZATION ////////////////////////////////////////////////// + //* Initialization of the Protection system (see CONTROL/Protection) + init_Prot_delay = 0.3; // Delay time initialization of the Protection devices (s) (All protections, expect the Chopper, are deactivated for the specified time ) + + /////////////////////////////////////////////////////////////_LOAD-FLOW CALCULATIONS ////////////////////////////////////////////////// + //* Load-Flow of the solar plant (see LF_init1) + LF_V = Vref_pu * Vgen_kVRMSLL * 1000 * sqrt(2 / 3); // Voltage for the LF device at the bus of the generator + + if (ControlMode == 1 || ControlMode == 3) { // if Q control or Q+Vac Control is selected, than PQ constraint is activated + LF_t1 = -1; // switch trick to active PQ constraint + LF_t2 = 1e15; // switch trick to active PQ constraint + LF_Ppark1_W = Spark_VA / (1.1162); // Active power in W of PQ constraint. The value is multiplied to account for losses and to remove 1.11 reserve factor also + LF_Ppark2_W = 0; // Active power in W of PV constraint. The value is set to zero + } else { // if Vac control is selected, than PV constraint is activated + LF_t1 = 1e15; // switch trick to active PV constraint + LF_t2 = -1; // switch trick to active PV constraint + LF_Ppark1_W = 0; // Active power in W of PQ constraint. The value is set to zero + LF_Ppark2_W = Spark_VA / (1.1162); // Active power in W of PV constraint. The value is multiplied to account for losses and to remove 1.11 reserve factor also + } + + + + + + // Q/V Regulator _ modified controller (updated version of VQ controller to calculate Iq_ref) + kp_var = 0.5; //0 + ki_var = 60; //40// + kp_volt = 1 // 0 + ki_volt = 60 // + + VoltReg_FRT_Gain = 2; // Proportional Gain of V control in FRT function + FRT_ON = 0.1; // FRT function is activated when voltage deviation from 1 pu (dV) is larger than FRT_ON + FRT_OFF = 0.075; // FRT function is deactivated when (dV < FRT_OFF) & (t_FRT > FRT_time) + FRT_time = 0.25; + + + + + + ////////////////////////////////////////////////////SOLAR PANEL PARAMETERS ////////////////////////////////// + /// ~75 W per module + // input + // Calcul + + //dVmpp=16.2; // maximum power voltage (V) under STC + //dnNbPVCellsSeries=Vdc_V / dVmpp; + //PoweratVmpp=73.9273; // new model diode + //dnNbPVCellsParallel=Pgen_W/dnNbPVCellsSeries/PoweratVmpp; + + Nmod_series = Vdc_V / VmaxP; // number of PV modules in series + Nmod_parallel = Ppark_W/(Vdc_V*ImaxP); // number of PV modules in parallel + //Nmod_parallel = Spark_VA/(Vdc_V*ImaxP); // number of PV modules in parallel + + ///////////// SOLAR PANEL PARAMETER CALCULATION from DATASHEET /////// (13/02/2015 Thomas Kauffmann) + + Vth = (Temp_ref + 273)*1.38065e-23/1.6022e-19; // diode threshold voltage + Ns = Ncell_series; + Io = Isc / (-1 + exp(Voc /(IdealFactor*Ns*Vth))); // diode reversed saturation current + Rs_max = Ns*((Voc - VmaxP)/(Ns*ImaxP) - exp(-Voc/(Vth*Ns*IdealFactor))*Vth*IdealFactor/Io); + + x1 = 10*Rs_max; // to initialize the while loop + x2 = Rs_max; + ite = 0; + + while ((x1 - x2)>num_tolerance){ + x1 = x2; + Rp_x1 = 1/(1/(VmaxP/ImaxP - x1) - Io/(IdealFactor*Ns*Vth)*exp((VmaxP + ImaxP*x1)/(IdealFactor*Ns*Vth))); + Iph_x1 = Isc*(VmaxP/(VmaxP - x1*ImaxP) - Io*x1*exp((VmaxP + ImaxP*x1)/(IdealFactor*Ns*Vth))/IdealFactor/Ns/Vth); + + // Newton Method + fx1 = Iph_x1 - ImaxP - (VmaxP + ImaxP*x1)/Rp_x1 - Io*(exp((VmaxP + ImaxP*x1)/(IdealFactor*Ns*Vth))-1); + df_x1 = ImaxP*VmaxP*(Isc - 2*ImaxP)/pow(VmaxP - x1*ImaxP,2) + Io*exp((VmaxP + ImaxP*x1)/(IdealFactor*Ns*Vth))*(ImaxP*VmaxP - Isc*IdealFactor*Ns*Vth + x1*ImaxP*(ImaxP - Isc))/pow(IdealFactor*Ns*Vth,2); + x2 = x1 - fx1/df_x1; + + ite = ite+1; + + } + + Rs = x2; + Rp = 1/(1/(VmaxP/ImaxP - x2) - Io/(IdealFactor*Ns*Vth)*exp((VmaxP + ImaxP*x2)/(IdealFactor*Ns*Vth))); + Iph = Isc*(Rs + Rp)/Rp; + + Rseries_PV = Rs*(Nmod_series/Nmod_parallel); + Rparallel_PV = Rp*(Nmod_series/Nmod_parallel); + + Vth_diode = (Temp + 273)*1.38065e-23/1.6022e-19; + I0_diode = Nmod_parallel*(Isc + Ki_pv*(Temp - Temp_ref))/(-1 + exp((Voc + Kv_pv*(Temp - Temp_ref))/(IdealFactor*Ns*Vth))); + Iph_T = Nmod_parallel*(Iph + Ki_pv*(Temp - Temp_ref)); + + Nseries = Ns*Nmod_series // total number of cell in series (value transmitted) + + +} diff --git a/EMTP-RV/ADN IPST2017/kth_IPST2017.ecf b/EMTP-RV/ADN IPST2017/kth_IPST2017.ecf new file mode 100644 index 0000000..724e3fe Binary files /dev/null and b/EMTP-RV/ADN IPST2017/kth_IPST2017.ecf differ