Simulation time increases considerably on implementing discrete PI controller

Question

I have developed a dynamic phasor model of a grid-connected single phase inverter. Running the open loop model of the simulation (with modulation index and initial angle constant) for a duration of 0.3 seconds with a step-size of 1e-05 takes about 45 seconds. However, implementing a discrete PI controller with sampling rate of 1e-4 leads to a simulation duration of about 30 minutes. My suspicion is the sample function in my PI code. Would be glad if anyone could help with identifying the cause and a possible workaround. The controller code is below. Thanks in advance.

model LinearController_DQ_InverterSP
 // reference currents
 parameter Real Id_ref[3] = {0.0, 10.0, 20.0}; 
 parameter Real Iq_ref[3] = {0.0, 5.0, 10.0};
 //step times
 parameter Real step_times[2] = {0.08, 0.12};
 // current reference currents
 Real Id_ref_curr;
 Real Iq_ref_curr;

 parameter Real Kp_cc = 0.0541 "Proportional term of controller";
 parameter Real Ki_cc = 55.577 "Integral term of controller";
 parameter Real T_sample = 1e-4 "Sampling period";
 
 // previous time step terms
 Real ud_prev(start = 0.0, fixed = true);
 Real uq_prev(start = 0.0, fixed = true);
 Real ed_prev(start = 0.0, fixed = true);
 Real eq_prev(start = 0.0, fixed = true);
 Real Vd_pcc_prev(start = 220*sqrt(2), fixed = true);
 Real Vq_pcc_prev(start = 0.0, fixed = true);
 
 // current time step terms
 Real ud_curr, uq_curr;
 Real ed_curr, eq_curr;
 Real Vin_d, Vin_q;

 Real a_const =  Kp_cc + Ki_cc * T_sample / 2;
 Real b_const = -Kp_cc + Ki_cc * T_sample / 2;
 
 Modelica.Blocks.Interfaces.RealOutput Mr(start = 0.87, fixed=true);  // modulation index
 Modelica.Blocks.Interfaces.RealOutput th(start=0.0, fixed=true) ; // angle
 Modelica.Blocks.Interfaces.RealInput Id ; // measured Id current 
 Modelica.Blocks.Interfaces.RealInput Iq ; // measure Iq current
 Modelica.Blocks.Interfaces.RealInput Vd_pcc ; // Vd of measured grid voltage
 Modelica.Blocks.Interfaces.RealInput Vq_pcc ; // Vq of measured grid voltage 


algorithm
   
   if time < step_times[1] then
     Id_ref_curr := Id_ref[1];
     Iq_ref_curr := Iq_ref[1];
   elseif (time > step_times[1]) and (time < step_times[2]) then
     Id_ref_curr := Id_ref[2];
     Iq_ref_curr := Iq_ref[2];
   else
     Id_ref_curr := Id_ref[3];
     Iq_ref_curr := Iq_ref[3];
   end if;
   
   when sample(0, T_sample) then
     //Calc error
       ed_curr := Id_ref_curr - Id;
       eq_curr := Iq_ref_curr - Iq;
     
    //Calc pid controller output
       ud_curr := ud_prev + a_const * ed_curr + b_const * ed_prev;
       uq_curr := uq_prev + a_const * eq_curr + b_const * eq_prev;
     
     //Apply feedforward term
       Vin_d := (Vd_pcc_prev + ud_prev);
       Vin_q := (Vq_pcc_prev + uq_prev);
     
     //Compute control signals
       th := atan(Vin_q / Vin_d);
       Mr := sqrt(Vin_d^2 + Vin_q^2) / 360;
     
     // Update values
       Vd_pcc_prev := Vd_pcc;
       Vq_pcc_prev := Vq_pcc;
       ed_prev := ed_curr;
       eq_prev := eq_curr;
       ud_prev := ud_curr;
       uq_prev := uq_curr;
   end when;

end LinearController_DQ_InverterSP;

Answer 1

Did you take a look at the step-size of the integrator and how it changes between the two scenarios? Usually, sampling (or any kind of discretization) forces the solver to take smaller steps by causing events. So e.g. enlarging the step-size from 1e-4s to 1e-3s should increase the performance by about a factor of 10.

Some references regarding this:

• https://mbe.modelica.university/behavior/discrete/events/
• http://people.inf.ethz.ch/~fcellier/Lect/NSDS/Ppt/nsds_21_engl.pdf.

Without having a closer look at the example (it doesn't work just by copy/paste), this would should be the explanation for your suspicion.

Answer 2

One issue is that you are trying to do the Modelica translator's job, by manually handling the previous variables, and writing the assignments in order in an algorithm.

The normal equation-oriented way would be:

model LinearController_DQ_InverterSP
 // reference currents
 parameter Real Id_ref[3] = {0.0, 10.0, 20.0};
 parameter Real Iq_ref[3] = {0.0, 5.0, 10.0};
 //step times
 parameter Real step_times[2] = {0.08, 0.12};
 // current reference currents
 Real Id_ref_curr;
 Real Iq_ref_curr;

 parameter Real Kp_cc = 0.0541 "Proportional term of controller";
 parameter Real Ki_cc = 55.577 "Integral term of controller";
 parameter Real T_sample = 1e-4 "Sampling period";


 // current time step terms
 Real ud_curr(start=0, fixed=true), uq_curr(start=0, fixed=true);
 Real ed_curr(start=0, fixed=true), eq_curr(start=0, fixed=true);
 Real Vin_d, Vin_q;

 Real a_const =  Kp_cc + Ki_cc * T_sample / 2;
 Real b_const = -Kp_cc + Ki_cc * T_sample / 2;

 Modelica.Blocks.Interfaces.RealOutput Mr(start = 0.87, fixed=true);  // modulation index
 Modelica.Blocks.Interfaces.RealOutput th(start=0.0, fixed=true);  // angle
 Modelica.Blocks.Interfaces.RealInput Id;  // measured Id current 
 Modelica.Blocks.Interfaces.RealInput Iq;  // measure Iq current
 Modelica.Blocks.Interfaces.RealInput Vd_pcc(start = 220*sqrt(2), fixed = true);  // Vd of measured grid voltage
 Modelica.Blocks.Interfaces.RealInput Vq_pcc(start = 0.0, fixed = true);  // Vq of measured grid voltage 


equation 

   if time < step_times[1] then
     Id_ref_curr =  Id_ref[1];
     Iq_ref_curr =  Iq_ref[1];
   elseif (time > step_times[1]) and (time < step_times[2]) then
     Id_ref_curr =  Id_ref[2];
     Iq_ref_curr =  Iq_ref[2];
   else
     Id_ref_curr =  Id_ref[3];
     Iq_ref_curr =  Iq_ref[3];
   end if;

   when sample(0, T_sample) then
     //Calc error
       ed_curr =  Id_ref_curr - Id;
       eq_curr =  Iq_ref_curr - Iq;

    //Calc pid controller output
       ud_curr =  pre(ud_curr) + a_const * ed_curr + b_const * pre(ed_curr);
       uq_curr =  pre(uq_curr) + a_const * eq_curr + b_const * pre(eq_curr);

     //Apply feedforward term
       Vin_d =  (pre(Vd_pcc) + pre(ud_curr));
       Vin_q =  (pre(Vq_pcc) + pre(uq_curr));

     //Compute control signals
       th =  atan(Vin_q / Vin_d);
       Mr =  sqrt(Vin_d^2 + Vin_q^2) / 360;
   end when;

end LinearController_DQ_InverterSP;

I haven't checked if there are odd interactions in terms of how the algorithm interacts with the assignments - but at least it is avoided in this variant.