Improvements + MMC model

- Initialization part in mask are now written in common scripts
- Tuning of Inertia Less control has been corrected
- Designation modification :
      GFo => GFM
      GFe => GFL
- MMC in GFM control has been added
- GFL input correction (dQ instead of dE)
- Power control loop is now a common block
- A tutorial has been added to explain how initialization is done

R2015b/Library/Init_Base.m

@@ -0,0 +1,18 @@
+%% Compute Base values
+% Can be used for all models
+% Index 1 = VSC side
+% Index 2 = grid side
+    fb = fn;
+    Pb = Pn;
+    cos_phi_n = Pn/Sn;
+    Sb = Pb/cos_phi_n;
+    Vn1 = Un1/sqrt(3);
+    Vb1 = Vn1;
+    Vn2 = Un2/sqrt(3);
+    Vb2 = Vn2;
+    Ib1 = Sb / (3*Vb1);
+    Ib2 = Sb / (3*Vb2);
+    Zb = (3*Vb1^2)/Sb;
+    wb = 2*pi*fb;  wb_pu=1;
+    Lb = Zb/wb;
+    Cb = 1/(Zb*wb);

R2015b/Library/Init_Cdc.m

@@ -0,0 +1,6 @@
+%% Initialization specific for model with an output DC Link
+
+% DC Link capacitor computation
+    C_DC = H_DC/1000*Pb/(1/2*Udcn^2);
+    set_param([gcb '/DC_Link_1'],'Capacitance',num2str(C_DC*2));
+    set_param([gcb '/DC_Link_2'],'Capacitance',num2str(C_DC*2));

R2015b/Library/Init_GFL.m

@@ -0,0 +1,11 @@
+%% Compute Control parameters for GFL VSC
+% Can be used for all GFL models
+
+% current control tuning
+    Lf_pu = Lt_pu;
+    Rf_pu = Rt_pu;
+    wn_i = 3/Tr_i;
+    Kp_i = 2*zeta_i*wn_i*Lf_pu/wb-Rf_pu;
+    Ti_i = 2*zeta_i/wn_i-Rf_pu/(wn_i^2*Lf_pu/wb);
+    Ki_i = Kp_i/Ti_i;
+    Kffi = 1;

R2015b/Library/Init_GFM.m

@@ -0,0 +1,38 @@
+%% Compute Control parameters for GFM VSC
+% Can be used for all GFM models
+
+% PLL derivation of the parameters
+    wn_PLL=5/Tr_PLL;
+    zeta_PLL = 1;
+    Kp_PLL_pu = 2*zeta_PLL*wn_PLL/wb;
+    Ki_PLL_pu = wn_PLL^2/wb;
+
+% Grid-Forming control tuning
+    % Set the chosen strategy
+        if (gainpopup == 1)
+            Chosen_Strategy = 1;
+        elseif (gainpopup == 2)
+            if (Param_Inertia == 1)
+                Chosen_Strategy = 2;
+            elseif (Param_Inertia == 2)
+                Chosen_Strategy = 3;
+            end
+        end
+    % Inertia Less
+    %wn_Inertia_Less = 2*pi*f_Inertia_Less;
+    %Ki_Inertia_Less = Lt_pu*wn_Inertia_Less^2/wb;
+    %Kp_Inertia_Less = 2*zeta_Inertia_Less*wn_Inertia_Less*Lt_pu/wb;
+    Kp_Inertia_Less = Lt_pu/wb / (Tr_Inertia_Less / 3000);
+
+    % With Inertia
+    Kc = 1/(Lt_pu);
+    kp_GFo = zeta_GFo*sqrt(2/(H_GFo*wb*Kc));
+    kd = 2*H_GFo*kp_GFo*Kc*wb;
+
+% Virtual Impedance tuning
+    X_VI = roots([1/DX_DR^2+1 2*Lt_pu Lt_pu^2-1/I_Max_VI_pu^2]);
+    X_VI = max(X_VI);
+    Kp_Rvi_pu = X_VI/(DX_DR*(I_Max_VI_pu-1)); 
+
+% Set the frequency droop gain
+    R = R_percent / 100;

R2015b/Library/Init_LF_L.m

@@ -0,0 +1,48 @@
+%% Initialization file for model with an output L-filter
+% Compute initial state of power and control component
+% Can be used for 2-level or MMC VSCs
+
+% Compute the delay compensation between control and power parts
+    if (Time_Step == -1)
+        Solver_Time_Step = 0;
+    else
+        Solver_Time_Step = Time_Step;
+    end
+
+% Compute Initial Voltage of VSC 
+%   from steady-state output transformer model and Loadflow result
+
+    R1_pu = Rt_pu/2;
+    L1_pu = Lt_pu/2;
+    R2_pu = Rt_pu/2;
+    L2_pu = Lt_pu/2;
+    X1_pu = L1_pu;
+    X2_pu = L2_pu;
+    Xm_pu = Lm_pu;
+
+    Zm_pu = 1/(1/Rm_pu+1/(1i*Xm_pu));
+
+    S0=P0+1i*Q0;
+    I2 = conj(S0)/(V0);
+
+    Vm = (R2_pu+1i*X2_pu)*I2 + V0;
+    Im = Vm/Zm_pu;
+
+    I1 = I2 + Im;
+    V1 = (R1_pu+1i*X1_pu) * I1 + Vm;
+
+    Vvsc0 = abs(V1);
+    Theta_vsc0 = angle(V1) * 180/pi+Theta0;
+
+    Va0 = Vvsc0*sqrt(2)*Vb1*sin(Theta_vsc0*pi/180);
+    Vb0 = Vvsc0*sqrt(2)*Vb1*sin(Theta_vsc0*pi/180-2*pi/3);
+    Vc0 = Vvsc0*sqrt(2)*Vb1*sin(Theta_vsc0*pi/180+2*pi/3);
+
+    Vm0_d = Vvsc0*cos(angle(V1));
+    Vm0_q = Vvsc0*sin(angle(V1));
+
+% Generate string to display LF results on Mask
+    mo = get_param(gcb,'MaskObject');
+    mo.getDialogControl('LF_vsc_str').Prompt = ['V0 = ', num2str(Vvsc0),' pu, Angle = ', num2str(Theta_vsc0),'°'];
+    mo.getDialogControl('LF_pcc_str').Prompt = ['V0 = ', num2str(V0),' pu, Angle = ', num2str(Theta0),'°'];
+    mo.getDialogControl('LF_pcc_Power_str').Prompt = ['P0 = ', num2str(P0),' pu, Q0 = ', num2str(Q0),' pu'];

R2015b/Library/Init_LF_LCL.m

@@ -0,0 +1,87 @@
+%% Initialization file for model with an output L-filter
+% Compute initial state of power and control component
+% Can be used for 2-level or MMC VSCs
+
+% Compute the delay compensation between control and power parts
+    if (Time_Step == -1)
+        Solver_Time_Step = 0;
+    else
+        Solver_Time_Step = Time_Step;
+    end
+
+% Compute Initial Voltage of VSC 
+%   from steady-state output transformer model and Loadflow result
+
+    % steady-state model of output transformer
+    R1_pu = Rt_pu/2;
+    L1_pu = Lt_pu/2;
+    R2_pu = Rt_pu/2;
+    L2_pu = Lt_pu/2;
+    X1_pu = L1_pu;
+    X2_pu = L2_pu;
+    Xm_pu = Lm_pu;
+    wb_pu = 1;
+
+    Zmag_pu = 1/(1/Rm_pu+1/(1i*Xm_pu));
+    S0=P0+1i*Q0;
+    I2 = conj(S0)/(V0);
+    Vmag = (R2_pu+1i*X2_pu)*I2 + V0;
+    Imag = Vmag/Zmag_pu;
+    I1 = I2 + Imag;
+    V1 = (R1_pu+1i*X1_pu) * I1 + Vmag;
+
+    Eg0 = abs(V1); % compute initial Eg value
+    Theta_m0 = angle(V1) * 180/pi+Theta0;
+
+    %Steady-state model of LC Filter
+    Is = I1 + 1i*Cf_pu*V1;
+    Theta_Is = angle(Is) * 180/pi+Theta0;
+    Vm = (1i*Lf_pu+Rf_pu)*Is+ V1;  %  
+
+% Generate string to display LF results on Mask
+
+    mo = get_param(gcb,'MaskObject');
+    mo.getDialogControl('LF_pcc_str').Prompt = ['V0 = ', num2str(V0),' pu, Angle = ', num2str(Theta0),'°'];
+    mo.getDialogControl('LF_pcc_Power_str').Prompt = ['P0 = ', num2str(P0),' pu, Q0 = ', num2str(Q0),' pu'];
+
+% Initial current and voltage states calculation 
+% Mainly for control loops initialization
+
+% Given Setpoint in PV bus
+
+    Vg0_d=V0*cos(Theta0*pi/180-Theta_m0*pi/180);
+    Vg0_q=V0*sin(Theta0*pi/180-Theta_m0*pi/180);
+    Vg0_dq = [Vg0_d ;Vg0_q; 0];
+
+% Initial state of AC grid current : Ig0dq
+
+    Theta_I2 = angle(I2) * 180/pi + Theta0;
+    Ig0_d = abs(I2)*cos(Theta_I2*pi/180-Theta_m0*pi/180);
+    Ig0_q = abs(I2)*sin(Theta_I2*pi/180-Theta_m0*pi/180);
+    Ig0_dq = [Ig0_d; Ig0_q; 0];
+
+% Initial state of capacitor voltage : eg0dq
+
+    Eg0_d = Eg0;
+    Eg0_q = 0;
+    Eg0_dq = [Eg0; 0; 0];
+
+% Initial state of VSC Current : Is0dq
+
+    Is0_d = abs(Is) * cos(Theta_Is*pi/180-Theta_m0*pi/180);
+    Is0_q = abs(Is) * sin(Theta_Is*pi/180-Theta_m0*pi/180);
+    Is0_dq = [Is0_d; Is0_q; 0];
+
+% Initial state of Modulation Voltage : Vm0dq
+
+    Vvsc0 = abs(Vm);
+    Theta_vsc0 = angle(Vm)*180/pi + Theta0;
+    Vm0_d = Vvsc0 * cos(Theta_vsc0*pi/180-Theta_m0*pi/180);
+    Vm0_q = Vvsc0 * sin(Theta_vsc0*pi/180-Theta_m0*pi/180);
+
+    Va0 = sqrt(2)*Vb1*Vvsc0*sin(Theta_vsc0*pi/180);
+    Vb0 = sqrt(2)*Vb1*Vvsc0*sin(Theta_vsc0*pi/180-2*pi/3);
+    Vc0 = sqrt(2)*Vb1*Vvsc0*sin(Theta_vsc0*pi/180+2*pi/3);
+    Vabc0 = [Va0 Vb0 Vc0];
+
+    mo.getDialogControl('LF_vsc_str').Prompt = ['V0 = ', num2str(Eg0),' pu, Angle = ', num2str(Theta_m0),'°'];

R2015b/Library/link_lib_sps.m

@@ -1,37 +1,86 @@
-model_list = libinfo(gcs,'IncludeCommented','on');
+function link_lib_sps(system_2_update)
+    model_list = libinfo(system_2_update,'IncludeCommented','on');
 
-find_lib = ["spsUniversalBridgeLib", 
-    "spsThreePhaseBreakerLib", 
-    "spsThreePhaseVIMeasurementLib", 
-    "spsSecondOrderFilterLib",
-    "spsSeriesRLCBranchLib",
-    "spsGroundLib",
-    "spsThreePhaseSeriesRLCBranchLib",
-    "spsSourcesModel",
-    "spsBreakerModel",
-    "spsThreePhaseTransformerTwoWindingsLib",
-    "spsThreePhaseSourceLib"];
+    find_lib = ["spsUniversalBridgeLib", 
+        "spsThreePhaseBreakerLib", 
+        "spsThreePhaseVIMeasurementLib", 
+        "spsSecondOrderFilterLib",
+        "spsSeriesRLCBranchLib",
+        "spsGroundLib",
+        "spsThreePhaseSeriesRLCBranchLib",
+        "spsSourcesModel",
+        "spsBreakerModel",
+        "spsThreePhaseTransformerTwoWindingsLib",
+        "spsThreePhaseSourceLib",
+        "spsDCVoltageSourceLib",
+        "spsLoadFlowBusLib",
+        "spsControlledVoltageSourceLib",
+        "spsBreakerLib",
+        "spsThreePhaseParallelRLCLoadLib",
+        "spspowerguiLib",
+        "spsVoltageMeasurementLib",
+        "spsPWMGenerator2LevelLib",
+        "spsThreePhaseFaultLib",
+        "spsThreePhaseDynamicLoadLib",
+        "spsAsynchronousMachinepuUnitsLib",
+        "spsThreePhasePISectionLineLib",
+        "spsThreePhaseSeriesRLCLoadLib",
+        "spsGroundingTransformerLib",
+        "spsDistributedParametersLineLib",
+        "spsSynchronousMachinepuStandardLib",
+        "spsHydraulicTurbineandGovernorLib",
+        "spsExcitationSystemLib",
+        "spsGenericPowerSystemStabilizerLib",
+        "spsMultiBandPowerSystemStabilizerLib",
+        "spsCurrentMeasurementLib",
+        "spsControlledCurrentSourceLib"];
 
-replace_lib = ["powerlib/Power Electronics",
-    "powerlib/Elements",
-    "powerlib/Measurements",
-    "powerlib_meascontrol/Filters",
-    "powerlib/Elements",
-    "powerlib/Elements",
-    "powerlib/Elements",
-    "powerlib/Electrical Sources",
-    "powerlib/Elements",
-    "powerlib/Elements",
-    "powerlib/Electrical Sources"];
+    replace_lib = ["powerlib/Power Electronics",
+        "powerlib/Elements",
+        "powerlib/Measurements",
+        "powerlib_meascontrol/Filters",
+        "powerlib/Elements",
+        "powerlib/Elements",
+        "powerlib/Elements",
+        "powerlib/Electrical Sources",
+        "powerlib/Elements",
+        "powerlib/Elements",
+        "powerlib/Electrical Sources",
+        "powerlib/Electrical Sources",
+        "powerlib/Measurements",
+        "powerlib/Electrical Sources",
+        "powerlib/Elements",
+        "powerlib/Elements",
+        "powerlib",
+        "powerlib/Measurements",
+        "powerlib_meascontrol/Pulse & Signal Generators",
+        "powerlib/Elements",
+        "powerlib/Elements",
+        "powerlib/Machines",
+        "powerlib/Elements",
+        "powerlib/Elements",
+        "powerlib/Elements",
+        "powerlib/Elements",
+        "powerlib/Machines",
+        "powerlib/Machines",
+        "powerlib/Machines",
+        "powerlib/Machines",
+        "powerlib/Machines",
+        "powerlib/Measurements",
+        "powerlib/Electrical Sources"];
 
-for i=1:length(model_list)
-    if contains(model_list(i).LinkStatus,'unresolved')
-        for j=1:length(find_lib)
-            if contains(model_list(i).Library,find_lib(j))
-                block_name = model_list(i).Block;
-                str = replace(string(model_list(i).ReferenceBlock),find_lib(j),replace_lib(j))
-                set_param(block_name,'SourceBlock',str);
-            end 
+    for i=1:length(model_list)
+        if contains(model_list(i).LinkStatus,'unresolved')
+            for j=1:length(find_lib)
+                if contains(model_list(i).Library,find_lib(j))
+                    block_name = model_list(i).Block;
+                    str = replace(string(model_list(i).ReferenceBlock),find_lib(j),replace_lib(j));
+                    if (j == 25)
+                        str = [str + ' '];
+                    end
+                    set_param(block_name,'SourceBlock',str);
+                end 
+            end
         end
     end
-end
+end

R2015b/Tutorial/01_Initialization/Init.m

@@ -0,0 +1,25 @@
+%% Grid Forming VSC with droop control strategies 
+
+Time_Step=10e-6; % Simulation time step (s) 
+
+fb = 50; % Nominal Frequency (Hz)
+Ug = 400e3; % Nominal Grid Voltage (V)
+
+%% Grid Parameters
+    wb = 2*pi*fb;
+
+% Nominal values  
+
+    Un = Ug;           % Phase-to-Phase nominal AC grid Voltage [V]  
+    Ub = Un;           % Base Phase-to-Phase Voltage value [V]  
+    Vb = Un/sqrt(3);   % Base simple Volatge value [V]    
+    Sb = 1e9;          % Base Apparent power [VA]  
+    Zb = Ub^2/Sb ;     % Base Grid Impedance  
+
+% Short circuit impedance
+
+    X_pu = 0.15;      % AC filter reactor in per-unit [p.u]
+    R_pu = 0.005;     % AC filter resistance in per-unit [p.u]
+    X    = X_pu*Zb;  % AC filter reactor in SI (International System of Units) 
+    L    = X/wb;     % AC filter inductance
+    R    = R_pu*Zb;  % AC filter resistance 

R2015b/Tutorial/01_Initialization/Solution/Solution.m

@@ -0,0 +1,8 @@
+% Inputs definition
+V1 = Un / sqrt(3);
+P = 500e6; Q = 100e6;
+
+% Quasi-static expression of V2
+% 3*V1*conj(I) = (P+jQ)
+
+V2 = (P-1i*Q)*(R+1i*X)/(3*V1)+V1;