Discrepancy in PI Controller Outputs for voltage control and current control for droop control strategy in MATLAB

Issue:

I am implementing a PI controller for the dq-transformation of a three-phase system in MATLAB. The goal is to control the q-axis voltage (iq_ref) and current (vq_inv) based on desired values. The system is simulated using a fixed-step solver in Simulink, and I am trying to replicate the results in MATLAB.

Problem Description:

I have implemented the PI controllers for Vq and Iq using the following equations:

% PI controller for Vq
error_vq(n) = Vq_ref(n) - V1q(n-1);
prop_vq = kpv * error_vq(n);
int_vq(n) = kiv * (Ts/2) * (error_vq(n) + error_vq(n-1)) + int_vq(n-1);
int_vq(n-1) = int_vq(n);
error_vq(n-1) = error_vq(n);
iq_ref(n) = prop_vq + int_vq(n) + (wc * V1d(n));

% PI controller for Iq
error_iq(n) = iq_ref(n) - I1q(n-1);
prop_iq = kpi * error_iq(n);
int_iq(n) = kii * (Ts/2) * (error_iq(n) + error_iq(n-1)) + int_iq(n-1);
int_iq(n-1) = int_iq(n);
error_iq(n-1) = error_iq(n);
Vq_inv(n) = prop_iq + int_iq(n) + (wL * I1d(n));

However, I noticed that the results obtained for iq_ref and Vq_inv in MATLAB are slightly different from the simulation results in Simulink. The values of these variables in MATLAB keep decreasing instead of remaining constant, as observed in the simulation.

Steps Taken:

1. Verified controller gains: The gains kpv, kiv, kpi, and kii have been double-checked and seem appropriate for the system.
2. Confirmed sampling time: The sampling time Ts used in the MATLAB code matches the actual sampling time used in the simulation.
3. Plotted variables: I plotted P_1, Q_1, and theta1 to understand their behavior during simulation, but I couldn't identify the issue.

Request:

I am seeking guidance on identifying the source of the discrepancy and how to modify the MATLAB code to match the simulation results from Simulink accurately. Any insights, suggestions, or potential fixes would be highly appreciated. Thank you in advance for your assistance.

Also, i can share .m and .slx file to understand this easily.

Images:

• whole simulation:

• from simulink to workspace:

• result from simulink of iq_ref:

• result from MATLAB program of iq_ref:

Code:

Va_1, Vb_1, Vc_1 and Ia_1, Ib_1, Ic_1 are imported from simulink in an array with sample time 1e-6 for 0.5 seconds

total samples = 500001

Ts=1e-06;
t = 0:Ts:0.5;
V1d=zeros(size(t));
V1q=zeros(size(t));
I1d=zeros(size(t));
I1q=zeros(size(t)); %  V1z=0.0;   I1z=0.0;
theta1=zeros(size(t));
Pf1new=zeros(size(t));  Pf1old=zeros(size(t));   Qf1new=zeros(size(t));   Qf1old=zeros(size(t));
P_1 = zeros(size(t));  Q_1 = zeros(size(t));

%f1error = [0.0 0.0];
%V1error; 
%V1out;
f1error = zeros(size(t));
error_vd= zeros(size(t));
error_vq= zeros(size(t));
error_id= zeros(size(t));
error_iq= zeros(size(t));
prop_vd=zeros(size(t));   prop_vq=zeros(size(t));   prop_id=zeros(size(t));   prop_iq=zeros(size(t));
int_vd = zeros(size(t));
int_vq = zeros(size(t));
int_id = zeros(size(t));
int_iq = zeros(size(t));
vanew=zeros(size(t));   vbnew=zeros(size(t));   vcnew=zeros(size(t));   vaold=zeros(size(t));   vbold=zeros(size(t));   vcold=zeros(size(t));
ianew=zeros(size(t));   ibnew=zeros(size(t));   icnew=zeros(size(t));   iaold=zeros(size(t));   ibold=zeros(size(t));   icold=zeros(size(t));
Vd_inv=zeros(size(t));   Vq_inv=zeros(size(t));   Vd_ref=zeros(size(t));   id_ref=zeros(size(t));   iq_ref=zeros(size(t));
m1=zeros(size(t));   n1=zeros(size(t));
V1d_off=zeros(size(t));   I1d_off=zeros(size(t));



wc=2*pi*50*9.94e-6; %0.0157; //314 rad/s*Cf
wL=2*pi*50*5.1e-3; %0.314; //314 rad/s*Lf

%//PI controller parameters
kpv=0.006245486195; %0.002
kiv=7.59150924004; %50
kpi=32.04424507; %2
kii=33740.7051; %50

%//Droop controller parameters
fnom = 50; %50.0
Vnom = 400; %10.0
Pbase = 10000*ones(size(t)); %500
%Qbase= 1*ones(size(t));
kf1=0.01; %0.01
kv1=0.05; %0.05

for n=2:500001
    vanew(n) = (Va_1(n)-vaold(1))+vaold(1);                  %  // Filter time constant is Ts*2pi*1000 for 1 kHz filter
    vaold(n) = vanew(n);
    
    vbnew(n) = (Vb_1(n)-vbold(1))+vbold(1);
    vbold(n) = vbnew(n);
    
    vcnew(n) = (Vc_1(n)-vcold(1))+vcold(1);
    vcold(n) = vcnew(n);
    
    ianew(n) = (Ia_1(n)-iaold(1))+iaold(1);
    iaold(n) = ianew(n);
    
    ibnew(n) = (Ib_1(n)-ibold(1))+ibold(1);
    ibold(n) = ibnew(n);
    
    icnew(n) = (Ic_1(n)-icold(1))+icold(1);
    icold(n) = icnew(n);
    
    V1d(n)=(2/3)*((sin(theta1(n)))*vanew(n)+(sin(theta1(n)-((2*pi)/3)))*vbnew(n)+(sin(theta1(n)+((2*pi)/3)))*vcnew(n));
    V1q(n)=(2/3)*((cos(theta1(n)))*vanew(n)+(cos(theta1(n)-((2*pi)/3)))*vbnew(n)+(cos(theta1(n)+((2*pi)/3)))*vcnew(n));
    I1d(n)=(2/3)*((sin(theta1(n)))*ianew(n)+(sin(theta1(n)-((2*pi)/3)))*ibnew(n)+(sin(theta1(n)+((2*pi)/3)))*icnew(n));
    I1q(n)=(2/3)*((cos(theta1(n)))*ianew(n)+(cos(theta1(n)-((2*pi)/3)))*ibnew(n)+(cos(theta1(n)+((2*pi)/3)))*icnew(n));
    
    
    
    P_1(n)=( V1d(n)*I1d(n) + V1q(n)*I1q(n) );
    Q_1(n)=( V1q(n)*I1d(n) - V1d(n)*I1q(n) );
    
    
    %/* Filter */
    Pf1new(n)=(P_1(n)-Pbase(n))+Pf1old(n);   %//Constant is Ts*wc=50e-6*31.41=0.0015705 0.0015705*
    Pf1old(n)=Pf1new(n);
    Qf1new(n)=(Q_1(n))+Qf1old(n);  % //Constant is Ts*wc=50e-6*31.41=0.0015705 0.0015705*
    Qf1old(n)=Qf1new(n);
    
    %//Droop control
    
    m1=kf1*(fnom/Pbase(n));
    n1=kv1*(Vnom/Pbase(n));
    
    %f1error(2) = (fnom-Pf1new(n)*m1) + f1error(1); %  0.0003142*         // Constant is 2*pi*Ts=2*3.14*50e-6=0.0003142
    %if(f1error(2)>6.28);
    %f1error(1) = f1error(2)-6.28; 
    %f1error(1) = f1error(2);
    %theta1(2)=f1error(1);
    %theta1(n+1)= (P_1(n)-Pbase)*theta1(1) + (2*pi*50*Ts); 
    theta1(n+1) = (2*pi*Ts)*(fnom-(Pf1new(n)*m1)) + theta1(n); %           // Constant is 2*pi*Ts=2*3.14*50e-6=0.0003142
    if theta1(n+1)>(2*pi)
        theta1(n+1) = theta1(n+1)-(2*pi);
        theta1(n+1) = theta1(n+1);
        theta1(n)=theta1(n+1);
    end
    
    Vd_ref(n) = (Vnom-(Qf1new(n)*n1));  %//Note that Vq_ref=0 by default
    Vq_ref(n)=0;
    % //PI control takes in Vd_ref and generates Vd_inv and Vq_inv
    
    % //PI controller for Vd
    error_vd(n)=((sqrt(2))*Vd_ref(n))-(V1d(n-1));
    prop_vd=kpv*error_vd(n);
    int_vd(n)=kiv*(Ts/2)*(error_vd(n)+error_vd(n-1))+int_vd(n-1);
    int_vd(n-1)=int_vd(n);
    error_vd(n-1)=error_vd(n);
    id_ref(n)=(prop_vd)+(int_vd(n))-(wc*V1q(n));         
    
    % //PI controller for Vq
    error_vq(n)=(Vq_ref(n))-(V1q(n-1));
    prop_vq=kpv*error_vq(n);
    int_vq(n)=kiv*(Ts/2)*(error_vq(n)+error_vq(n-1))+int_vq(n-1);
    int_vq(n-1)=int_vq(n);
    error_vq(n-1)=error_vq(n);
    iq_ref(n)=(prop_vq)+(int_vq(n))+(wc*V1d(n));
    
    %//PI controller for Id
    error_id(n)=id_ref(n)-(I1d(n-1));
    prop_id=kpi*error_id(n);
    int_id(n)=kii*(Ts/2)*(error_id(n)+error_id(n-1))+int_id(n-1);
    int_id(n-1)=int_id(n);
    error_id(n-1)=error_id(n);
    Vd_inv(n)=prop_id+int_id(n)-(wL*I1q(n));
    
   %//PI controller for Iq
    error_iq(n)=iq_ref(n)-(I1q(n-1));
    prop_iq=kpi*error_iq(n);
    int_iq(n)=kii*(Ts/2)*(error_iq(n)+error_iq(n-1))+int_iq(n-1);
    int_iq(n-1)=int_iq(n);
    error_iq(n-1)=error_iq(n);
    Vq_inv(n)=prop_iq+int_iq(n)+(wL*I1d(n));
    
    
end

[Sam Chak]

Hi @Abdul Moiz Qureshi

I am unable to plot iq_ref for comparison purposes because the input values, Va_1, Vb_1, Vc_1 and Ia_1, Ib_1, Ic_1 (from Simulink), are missing. You need to provide the data, maybe in mat-file.

Ts=1e-06;
t = 0:Ts:0.5;
V1d=zeros(size(t));
V1q=zeros(size(t));
I1d=zeros(size(t));
I1q=zeros(size(t)); %  V1z=0.0;   I1z=0.0;
theta1=zeros(size(t));
Pf1new=zeros(size(t));  Pf1old=zeros(size(t));   Qf1new=zeros(size(t));   Qf1old=zeros(size(t));
P_1 = zeros(size(t));  Q_1 = zeros(size(t));
%f1error = [0.0 0.0];
%V1error; 
%V1out;
f1error = zeros(size(t));
error_vd= zeros(size(t));
error_vq= zeros(size(t));
error_id= zeros(size(t));
error_iq= zeros(size(t));
prop_vd=zeros(size(t));   prop_vq=zeros(size(t));   prop_id=zeros(size(t));   prop_iq=zeros(size(t));
int_vd = zeros(size(t));
int_vq = zeros(size(t));
int_id = zeros(size(t));
int_iq = zeros(size(t));
vanew=zeros(size(t));   vbnew=zeros(size(t));   vcnew=zeros(size(t));   vaold=zeros(size(t));   vbold=zeros(size(t));   vcold=zeros(size(t));
ianew=zeros(size(t));   ibnew=zeros(size(t));   icnew=zeros(size(t));   iaold=zeros(size(t));   ibold=zeros(size(t));   icold=zeros(size(t));
Vd_inv=zeros(size(t));   Vq_inv=zeros(size(t));   Vd_ref=zeros(size(t));   id_ref=zeros(size(t));   iq_ref=zeros(size(t));
m1=zeros(size(t));   n1=zeros(size(t));
V1d_off=zeros(size(t));   I1d_off=zeros(size(t));
wc=2*pi*50*9.94e-6; %0.0157; //314 rad/s*Cf
wL=2*pi*50*5.1e-3; %0.314; //314 rad/s*Lf
%//PI controller parameters
kpv=0.006245486195; %0.002
kiv=7.59150924004; %50
kpi=32.04424507; %2
kii=33740.7051; %50
%//Droop controller parameters
fnom = 50; %50.0
Vnom = 400; %10.0
Pbase = 10000*ones(size(t)); %500
%Qbase= 1*ones(size(t));
kf1=0.01; %0.01
kv1=0.05; %0.05
for n=2:500001
    vanew(n) = (Va_1(n)-vaold(1))+vaold(1);                  %  // Filter time constant is Ts*2pi*1000 for 1 kHz filter
    vaold(n) = vanew(n);
    
    vbnew(n) = (Vb_1(n)-vbold(1))+vbold(1);
    vbold(n) = vbnew(n);
    
    vcnew(n) = (Vc_1(n)-vcold(1))+vcold(1);
    vcold(n) = vcnew(n);
    
    ianew(n) = (Ia_1(n)-iaold(1))+iaold(1);
    iaold(n) = ianew(n);
    
    ibnew(n) = (Ib_1(n)-ibold(1))+ibold(1);
    ibold(n) = ibnew(n);
    
    icnew(n) = (Ic_1(n)-icold(1))+icold(1);
    icold(n) = icnew(n);
    
    V1d(n)=(2/3)*((sin(theta1(n)))*vanew(n)+(sin(theta1(n)-((2*pi)/3)))*vbnew(n)+(sin(theta1(n)+((2*pi)/3)))*vcnew(n));
    V1q(n)=(2/3)*((cos(theta1(n)))*vanew(n)+(cos(theta1(n)-((2*pi)/3)))*vbnew(n)+(cos(theta1(n)+((2*pi)/3)))*vcnew(n));
    I1d(n)=(2/3)*((sin(theta1(n)))*ianew(n)+(sin(theta1(n)-((2*pi)/3)))*ibnew(n)+(sin(theta1(n)+((2*pi)/3)))*icnew(n));
    I1q(n)=(2/3)*((cos(theta1(n)))*ianew(n)+(cos(theta1(n)-((2*pi)/3)))*ibnew(n)+(cos(theta1(n)+((2*pi)/3)))*icnew(n));
    
    
    
    P_1(n)=( V1d(n)*I1d(n) + V1q(n)*I1q(n) );
    Q_1(n)=( V1q(n)*I1d(n) - V1d(n)*I1q(n) );
    
    
    %/* Filter */
    Pf1new(n)=(P_1(n)-Pbase(n))+Pf1old(n);   %//Constant is Ts*wc=50e-6*31.41=0.0015705 0.0015705*
    Pf1old(n)=Pf1new(n);
    Qf1new(n)=(Q_1(n))+Qf1old(n);  % //Constant is Ts*wc=50e-6*31.41=0.0015705 0.0015705*
    Qf1old(n)=Qf1new(n);
    
    %//Droop control
    
    m1=kf1*(fnom/Pbase(n));
    n1=kv1*(Vnom/Pbase(n));
    
    %f1error(2) = (fnom-Pf1new(n)*m1) + f1error(1); %  0.0003142*         // Constant is 2*pi*Ts=2*3.14*50e-6=0.0003142
    %if(f1error(2)>6.28);
    %f1error(1) = f1error(2)-6.28; 
    %f1error(1) = f1error(2);
    %theta1(2)=f1error(1);
    %theta1(n+1)= (P_1(n)-Pbase)*theta1(1) + (2*pi*50*Ts); 
    theta1(n+1) = (2*pi*Ts)*(fnom-(Pf1new(n)*m1)) + theta1(n); %           // Constant is 2*pi*Ts=2*3.14*50e-6=0.0003142
    if theta1(n+1)>(2*pi)
        theta1(n+1) = theta1(n+1)-(2*pi);
        theta1(n+1) = theta1(n+1);
        theta1(n)=theta1(n+1);
    end
    
    Vd_ref(n) = (Vnom-(Qf1new(n)*n1));  %//Note that Vq_ref=0 by default
    Vq_ref(n)=0;
    % //PI control takes in Vd_ref and generates Vd_inv and Vq_inv
    
    % //PI controller for Vd
    error_vd(n)=((sqrt(2))*Vd_ref(n))-(V1d(n-1));
    prop_vd=kpv*error_vd(n);
    int_vd(n)=kiv*(Ts/2)*(error_vd(n)+error_vd(n-1))+int_vd(n-1);
    int_vd(n-1)=int_vd(n);
    error_vd(n-1)=error_vd(n);
    id_ref(n)=(prop_vd)+(int_vd(n))-(wc*V1q(n));         
    
    % //PI controller for Vq
    error_vq(n)=(Vq_ref(n))-(V1q(n-1));
    prop_vq=kpv*error_vq(n);
    int_vq(n)=kiv*(Ts/2)*(error_vq(n)+error_vq(n-1))+int_vq(n-1);
    int_vq(n-1)=int_vq(n);
    error_vq(n-1)=error_vq(n);
    iq_ref(n)=(prop_vq)+(int_vq(n))+(wc*V1d(n));
    
    %//PI controller for Id
    error_id(n)=id_ref(n)-(I1d(n-1));
    prop_id=kpi*error_id(n);
    int_id(n)=kii*(Ts/2)*(error_id(n)+error_id(n-1))+int_id(n-1);
    int_id(n-1)=int_id(n);
    error_id(n-1)=error_id(n);
    Vd_inv(n)=prop_id+int_id(n)-(wL*I1q(n));
    
    % //PI controller for Vq
    error_vq(n)=(Vq_ref(n))-(V1q(n-1));
    prop_vq=kpv*error_vq(n);
    int_vq(n)=kiv*(Ts/2)*(error_vq(n)+error_vq(n-1))+int_vq(n-1);
    int_vq(n-1)=int_vq(n);
    error_vq(n-1)=error_vq(n);
    iq_ref(n)=(prop_vq)+(int_vq(n))+(wc*V1d(n));
    
    
end

Unrecognized function or variable 'Va_1'.

[Abdul Moiz Qureshi]

variables.mat

[Sam Chak]

Hi @Abdul Moiz Qureshi

Thanks for the data. From the equation for iq_ref, (wc = 0.0031 is a relatively small constant)

iq_ref(n) = kpv*error_vq(n) + int_vq(n) + wc*V1d(n);

and the tiled plots below, we can see that int_vq is responsible for shaping how iq_ref behaves. Thus, it is worth checking the logical execution of these two lines:

int_vq(n) = kiv*(Ts/2)*(error_vq(n) + error_vq(n-1)) + int_vq(n-1);

int_vq(n-1) = int_vq(n);

Is the term int_vq(n-1) causing the decrement in int_vq?

[Abdul Moiz Qureshi]

yes I just checked that yesterday, int_vq(n-1) is responsible for decrement but it is responsible for the amplitude same as the simulation too. Is there any way to fix that? like getting the same plot without decrement.

Thanks :)

Answer 1

[Paul]

With regard to this part of the code:

    
    % //PI controller for Vd
    error_vd(n)=((sqrt(2))*Vd_ref(n))-(V1d(n-1));
    prop_vd=kpv*error_vd(n);
    int_vd(n)=kiv*(Ts/2)*(error_vd(n)+error_vd(n-1))+int_vd(n-1);
    int_vd(n-1)=int_vd(n);
    error_vd(n-1)=error_vd(n);
    id_ref(n)=(prop_vd)+(int_vd(n))-(wc*V1q(n));         
    
    % //PI controller for Vq
    error_vq(n)=(Vq_ref(n))-(V1q(n-1));
    prop_vq=kpv*error_vq(n);
    int_vq(n)=kiv*(Ts/2)*(error_vq(n)+error_vq(n-1))+int_vq(n-1);
    int_vq(n-1)=int_vq(n);
    error_vq(n-1)=error_vq(n);
    iq_ref(n)=(prop_vq)+(int_vq(n))+(wc*V1d(n));
    
    %//PI controller for Id
    error_id(n)=id_ref(n)-(I1d(n-1));
    prop_id=kpi*error_id(n);
    int_id(n)=kii*(Ts/2)*(error_id(n)+error_id(n-1))+int_id(n-1);
    int_id(n-1)=int_id(n);
    error_id(n-1)=error_id(n);
    Vd_inv(n)=prop_id+int_id(n)-(wL*I1q(n));
    
    % //PI controller for Vq
    error_vq(n)=(Vq_ref(n))-(V1q(n-1));
    prop_vq=kpv*error_vq(n);
    int_vq(n)=kiv*(Ts/2)*(error_vq(n)+error_vq(n-1))+int_vq(n-1);
    int_vq(n-1)=int_vq(n);
    error_vq(n-1)=error_vq(n);
    iq_ref(n)=(prop_vq)+(int_vq(n))+(wc*V1d(n));
    

Why is the section for "PI controller for Vq" implemented twice?

What is the purpose of these lines in the PI controller for Vd section

int_vd(n-1)=int_vd(n);
error_vd(n-1)=error_vd(n);

and similarly for the other sections?

[Abdul Moiz Qureshi]

The purpose of the lines in the PI controller for Vd section, which are:

• int_vd(n-1)=int_vd(n);
• error_vd(n-1)=error_vd(n);
• Integral (I) Term Update: The line int_vd(n-1)=int_vd(n); updates the integral component (I-term) of the PI controller for Vd. It assigns the current value of the integral term int_vd(n) to the previous time step value int_vd(n-1). This update ensures that the previous integral state is maintained and carried forward for the next iteration.
• Error Update: The line error_vd(n-1)=error_vd(n); updates the error value used in the next iteration of the controller. It assigns the current error value error_vd(n) to the previous time step value error_vd(n-1). This update is essential to have the correct error value for the integral calculation in the next time step.

In summary, these lines are responsible for updating the integral term and the error value of the PI controller at each time step, ensuring proper functioning and accumulation of integral error for effective control of the system.

also I mistakenly wrote twice Vq, actually it is for 'Iq' and now i updated it you can check.

[Paul]

Except that the variables int_vd and error_vd are being stored as arrays as each element is computed, so accessing the previous value in the current iteration should be accomplished with indexing. It actually looks like those assignments don't do anyhting, unless those variables are used downstream for some some other function.

Consider the situation when n = 3. int_vd(2) and error_vd(2) are used to update int_vd(3), and then they are reassigned.

int_vd(3)=kiv*(Ts/2)*(error_vd(3)+error_vd(2))+int_vd(2);
int_vd(2)=int_vd(3);
error_vd(2)=error_vd(3);

Now, on the next pass n = 4

int_vd(4)=kiv*(Ts/2)*(error_vd(4)+error_vd(3))+int_vd(3);
int_vd(3)=int_vd(4);
error_vd(3)=error_vd(4);

int_vd(4) is updated based on int_vd(3) and error_vd(3), which means the values assigned to int_vd(2) and error_vd(2) on the n=3 iteration do not mean anything relative to computing the integral.

Don't know if this is causing any problems in the code, but maybe is an indication of incorrect logic.

Also, I'm not sure that the calculations of the integrals aren't off by one time step. Would have to look at that in a bit more detail.

[Abdul Moiz Qureshi]

how it can be updated without 'n' mention in the bracket for iterations like you share the particular number of iteration for understanding or should i mention n+1 or n?

also the situation is that I've to measure theta1 first so I've to calculate it for '0' first cause as theta1 is important for dq-transformation and also for power calculations but first dq-transformation then power calculations as they both are synchronise or utilizing theta1 at the same time.

And if you want to know more I can share with you.

Thanks

[Paul]

If you're going to store all of the relevant variables in arrays, then of course you need to index with n in the equations internal to the iteration. I was only trying to show that the assignments to the (n-1) variables on each step through the iteration seem to be not necessary for computing the PI control.

[Abdul Moiz Qureshi]

I understand that still getting a different value. Thanks for correction but unfortunately it can't help as such. :(

[Abdul Moiz Qureshi]

@Paul can you please help me out?

[Paul]

Another thing I noticed ....

In these lines

 % //PI controller for Vd
    error_vd(n)=((sqrt(2))*Vd_ref(n))-(V1d(n-1));
    prop_vd=kpv*error_vd(n);
    int_vd(n)=kiv*(Ts/2)*(error_vd(n)+error_vd(n-1))+int_vd(n-1);
    int_vd(n-1)=int_vd(n);
    error_vd(n-1)=error_vd(n);
    id_ref(n)=(prop_vd)+(int_vd(n))-(wc*V1q(n));         
    
    % //PI controller for Vq
    error_vq(n)=(Vq_ref(n))-(V1q(n-1));
    prop_vq=kpv*error_vq(n);
    int_vq(n)=kiv*(Ts/2)*(error_vq(n)+error_vq(n-1))+int_vq(n-1);
    int_vq(n-1)=int_vq(n);
    error_vq(n-1)=error_vq(n);
    iq_ref(n)=(prop_vq)+(int_vq(n))+(wc*V1d(n));
    
    %//PI controller for Id
    error_id(n)=id_ref(n)-(I1d(n-1));
    prop_id=kpi*error_id(n);
    int_id(n)=kii*(Ts/2)*(error_id(n)+error_id(n-1))+int_id(n-1);
    int_id(n-1)=int_id(n);
    error_id(n-1)=error_id(n);
    Vd_inv(n)=prop_id+int_id(n)-(wL*I1q(n));
    
   %//PI controller for Iq
    error_iq(n)=iq_ref(n)-(I1q(n-1));
    prop_iq=kpi*error_iq(n);
    int_iq(n)=kii*(Ts/2)*(error_iq(n)+error_iq(n-1))+int_iq(n-1);
    int_iq(n-1)=int_iq(n);
    error_iq(n-1)=error_iq(n);
    Vq_inv(n)=prop_iq+int_iq(n)+(wL*I1d(n));
    

In all four sections, the error signal is formed by: var(n) - anothervar(n-1). The different indices looks suspicious. Also, idref(n) and iqref(n) are computed in the first two sections and then each are used to compute an error signal in the following sections. That also looks a bit odd.

Can't say anything definitively though. Even after zoomin in on the block diagram, I wasn't sure which part of the diagram these equations are supposed to implement.

[Abdul Moiz Qureshi]

The equations are related to voltage and current control block. I can share the .slx file too.

[Abdul Moiz Qureshi]

Hey @Paul, this is what i trying to implement hope this can help you to understand the problem in detail

The simulation is based on three phase inverter with droop control strategy.

Where three phase load voltages and current first converted into dq-transformation Vd, Vq, Id, Iq but I don't have theta initially as theta is found by the active power as in droop control P-f so initially compute both as theta=zero after that theta is calculated it is continue to increment for the no. of samples (500001 ) imported from simulink to MATLAB workspace. The Vd_ref is also computed from Q-v relation in droop control.

After that there is voltage control where Vd_ref is compare with Vd (from load) then passes through PI controller then Vd*(2*pi*f*c) subtracted to id_ref after that same goes for Vq_ref (0) to compute iq_ref.

Then id_ref and iq_ref passes in current control and compare with Id and Iq (from load) respectively then into PI controller further into Id*(2*pi*f*L) subtracted from it to get Vd_inv and same for iq_ref to get Vq_inv. And Vd_inv, Vq_inv both passes through inverse dq-transform here also theta use to get three phase voltages references to generate SPWM.

Remember that the first iteration as there is no theta available for transformation then it is taken as zero then vd,vq, id, iq used to calculate active power which is responsible to generate theta then this theta on next iteration used to convert three phase voltages and currents into dq-transform and then power to theta and the cycle goes on till the no. of samples.

I just import the Va, Vb, Vc and Ia, Ib, Ic which are responsible to transform into dq-axis and then into active and reactive but for both I need theta which is comes from P-f droop control and before that I don't have any theta so at first take it as zero and then after conversion into dq and then power and from power we can get frequency and theta so that theta is utilized for next iteration.Where three phase load voltages and current first converted into dq-transformation Vd, Vq, Id, Iq but I don't have theta initially as theta is found by the active power as in droop control P-f so initially compute both as theta=zero after that theta is calculated it is continue to increment for the no. of samples (500001 ) imported from simulink to MATLAB workspace. The Vd_ref is also computed from Q-v relation in droop control.After that there is voltage control where Vd_ref is compare with Vd (from load) then passes through PI controller then Vd*(2*pi*f*c) subtracted to id_ref after that same goes for Vq_ref (0) to compute iq_ref.Then id_ref and iq_ref passes in current control and compare with Id and Iq (from load) respectively then into PI controller further into Id*(2*pi*f*L) subtracted from it to get Vd_inv and same for iq_ref to get Vq_inv. And Vd_inv, Vq_inv both passes through inverse dq-transform here also theta use to get three phase voltages references to generate SPWM. Remember that the first iteration as there is no theta available for transformation then it is taken as zero then vd,vq, id, iq used to calculate active power which is responsible to generate theta then this theta on next iteration used to convert three phase voltages and currents into dq-transform and then power to theta and the cycle goes on till the no. of samples.I just import the Va, Vb, Vc and Ia, Ib, Ic which are responsible to transform into dq-axis and then into active and reactive but for both I need theta which is comes from P-f droop control and before that I don't have any theta so at first take it as zero and then after conversion into dq and then power and from power we can get frequency and theta so that theta is utilized for next iteration.

Answer 2

[Sam Chak]

Hi @Abdul Moiz Qureshi

I do not understand the logical operations of the two lines, because they are in code form and my mind is not good at translating code into mathematical equations. But when I adjust the coefficient of int_vq(n-1) from 1 to 0.9999995, it affects iq_ref.

int_vq(n)   = kiv*(Ts/2)*(error_vq(n) + error_vq(n-1)) + 0.9999995*int_vq(n-1);
int_vq(n-1) = int_vq(n);

I suggest you consult with @Paul if he has any idea on the math of int_vq(n).

[Abdul Moiz Qureshi]

As you know that int_vq(n-1) responsible for the dip in a plot but why it works fine for id_ref but not for iq_ref?

[Abdul Moiz Qureshi]

@Sam Chak can you please check the anti-windup tweaking i did here!

Ts=1e-06;
t = 0:Ts:0.5;

theta1=[0.0,0.0];
P_1 = 0.0;  Q_1 = 0.0;
P_2 = 0; Q_2 =  0;
frequency = [0,0];

  prop_vd = 0.0;   prop_vq = 0.0;
  prop_id = 0.0;   prop_iq = 0.0;
    
  Vd_inv = 0.0;  Vq_inv= 0.0;  
  Vd_ref = 0.0;  Vq_ref= 0.0;
  id_ref = 0.0;  iq_ref= 0.0;
 
  sine_a= 0; sine_b= 0; sine_c = 0;

  V1d= 0.0; V1q= 0.0; I1d= 0.0; I1q= 0.0;
 
  error_vd =  [0.0,0.0];
  error_vq =  [0.0,0.0];
  error_id =  [0.0,0.0];
  error_iq =  [0.0,0.0];
  
  int_vd =  [0.0,0.0];
  int_vq =  [0.0,0.0];
  int_id =  [0.0,0.0];
  int_iq =  [0.0,0.0];

wc= 3.122743098e-3;  %0.0157; //314 rad/s*Cf %0.003122743098;
wL= 1.602212253e0; %0.314; //314 rad/s*Lf %1.602212253;

%//PI controller parameters
kpv=0.006245486195; %0.002
kiv=7.59150924004; %50
kpi=32.04424507; %2
kii=33740.7051; %50

%//Droop controller parameters
fnom = 50; %50.0
Vnom = 400; %10.0
Pbase = 10000; %500

kf1=0.01; %0.01
kv1=0.05; %0.05


for n=2:500001
    
    %dt = t(n) - t(n-1);

    V1d(n)=(2/3)*((sin(theta1(n)))*Va_1(n)+(sin(theta1(n)-((2*pi)/3)))*Vb_1(n)+(sin(theta1(n)+((2*pi)/3)))*Vc_1(n));
    V1q(n)=(2/3)*((cos(theta1(n)))*Va_1(n)+(cos(theta1(n)-((2*pi)/3)))*Vb_1(n)+(cos(theta1(n)+((2*pi)/3)))*Vc_1(n));
    I1d(n)=(2/3)*((sin(theta1(n)))*Ia_1(n)+(sin(theta1(n)-((2*pi)/3)))*Ib_1(n)+(sin(theta1(n)+((2*pi)/3)))*Ic_1(n));
    I1q(n)=(2/3)*((cos(theta1(n)))*Ia_1(n)+(cos(theta1(n)-((2*pi)/3)))*Ib_1(n)+(cos(theta1(n)+((2*pi)/3)))*Ic_1(n));



   P_1(n)= V1d(n)*I1d(n) + V1q(n)*I1q(n) ;
   Q_1(n)= V1q(n)*I1d(n) - V1d(n)*I1q(n) ;


 %/*  lter */
    P_2(n)=(P_1(n)-Pbase);   %//Constant is Ts*wc=50e-6*31.41=0.0015705 0.0015705*
    %Pf1old(n)=P_2;
    Q_2(n)=(Q_1(n));  % //Constant is Ts*wc=50e-6*31.41=0.0015705 0.0015705*
    %Qf1old=Q_2;
    
 %//Droop control

   m1=kf1*(fnom/Pbase);
   n1=kv1*(Vnom/Pbase);

  
%theta1(n+1) = (2*pi*Ts)*(f1error-2)
    

   Vd_ref(n) = (Vnom-(Q_2(n)*n1));  %//Note that Vq_ref=0 by default
   Vq_ref(n)=0;

   % //PI control takes in Vd_ref and generates Vd_inv and Vq_inv

   % //PI controller for Vd
    error_vd(n)=sqrt(2) * Vd_ref(n) - V1d(n);
    prop_vd(n)=kpv * error_vd(n);
    int_vd(n)= int_vd(n-1) + kiv * Ts * error_vd(n);
    if int_vd(n) > 1
        int_vd(n+1) = 1;
    elseif int_vd(n) < -1
        int_vd(n+1) = -1;
    end
    id_ref(n)=prop_vd(n) + int_vd(n) - wc * V1q(n);         
    
    % //PI controller for Vq
    error_vq(n)= -V1q(n);
    prop_vq(n)=kpv * error_vq(n);
    int_vq(n)=int_vq(n-1) + kiv * Ts * error_vq(n);
    if int_vq(n) > 1
        int_vq(n+1) = 1;
    elseif int_vq(n) < -1
        int_vq(n+1) = -1;
    end
    iq_ref(n)=prop_vq(n) + int_vq(n) + wc * V1d(n);
    
    %//PI controller for Id
    error_id(n)=id_ref(n) - I1d(n);
    prop_id(n)=kpi * error_id(n);
    int_id(n)=int_id(n-1) + kii * Ts * error_id(n);
    if int_id(n) > 1
        int_id(n+1) = 1;
    elseif int_id(n) < -1
        int_id(n+1) = -1;
    end
    Vd_inv(n)=prop_id(n) + int_id(n) - wL * I1q(n);
    
   %//PI controller for Iq
    error_iq(n)=iq_ref(n) - I1q(n);
    prop_iq(n)=kpi * error_iq(n);
    int_iq(n)=int_iq(n-1) + kii* Ts * error_iq(n);
    if int_iq(n) > 1
        int_iq(n+1) = 1;
    elseif int_iq(n) < -1
        int_iq(n+1) = -1;
    end
    Vq_inv(n)=prop_iq(n) + int_iq(n) + wL*I1d(n);
    
   frequency = (fnom-(P_2(n)*m1)); %  0.0003142*         // Constant is 2*pi*Ts=2*3.14*50e-6=0.0003142
   %theta1(n+1) = (2*pi*Ts)*(fnom-(P_2(n)*m1)) + theta1(n); %           // Constant is 2*pi*Ts=2*3.14*50e-6=0.0003142
   theta1(n+1)=theta1(n)+(2*pi*(frequency*Ts));
%if theta1(n+1)>=((2*pi)-(2*pi)*(frequency*Ts))
%theta1(n+1)=0;
%end
 %//Inverse dq transform
   sine_a(n)=Vd_inv(n)*sin(theta1(n))+Vq_inv(n)*cos(theta1(n));
   sine_b(n)=Vd_inv(n)*sin(theta1(n))-((2*pi)/3)+Vq_inv(n)*cos(theta1(n)-((2*pi)/3));
   sine_c(n)=Vd_inv(n)*sin(theta1(n))+((2*pi)/3)+Vq_inv(n)*cos(theta1(n)+((2*pi)/3));


end

   figure;
    subplot(2, 2, 1);
    plot(V1d);
    title('V1d');
    subplot(2, 2, 2);
    plot(V1q);
    title('V1q');
    subplot(2, 2, 3);
    plot(I1d);
    title('I1d');
    subplot(2, 2, 4);
    plot(I1q);
    title('I1q');

     figure;
    plot(theta1);
    title('theta wt');

     figure;
    subplot(2, 1, 1);
    plot(P_1);
    title('Active Power P_1');
    subplot(2, 1, 2);
    plot(Q_1);
    title('Reactive Power Q_1');

     figure;
    subplot(2, 1, 1);
    plot(Vd_inv);
    title('Vd_inv');
    subplot(2, 1, 2);
    plot(Vq_inv);
    title('Vq_inv');

  figure;
    subplot(2, 1, 1);
    plot(id_ref);
    title('id_ref');
    subplot(2, 1, 2);
    plot(iq_ref);
    title('iq_ref');

     figure;
    plot(sine_a)
    title('Inverse DQ')
    hold on
    plot(sine_b)
    hold on
    plot(sine_c)

 figure;   
plot(V1d)
   hold on
    plot(V1q)
    grid ont = 0:Ts:0.5;theta1=[0.0,0.0];P_1 = 0.0;  Q_1 = 0.0;P_2 = 0; Q_2 =  0;frequency = [0,0];prop_vd = 0.0;   prop_vq = 0.0;prop_id = 0.0;   prop_iq = 0.0;  Vd_inv = 0.0;  Vq_inv= 0.0;  Vd_ref = 0.0;  Vq_ref= 0.0;id_ref = 0.0;  iq_ref= 0.0; sine_a= 0; sine_b= 0; sine_c = 0;V1d= 0.0; V1q= 0.0; I1d= 0.0; I1q= 0.0; error_vd =  [0.0,0.0];error_vq =  [0.0,0.0];error_id =  [0.0,0.0];error_iq =  [0.0,0.0];int_vd =  [0.0,0.0];int_vq =  [0.0,0.0];int_id =  [0.0,0.0];int_iq =  [0.0,0.0];wc= 3.122743098e-3;  %0.0157; //314 rad/s*Cf %0.003122743098;wL= 1.602212253e0; %0.314; //314 rad/s*Lf %1.602212253;%//PI controller parameterskpv=0.006245486195; %0.002kiv=7.59150924004; %50kpi=32.04424507; %2kii=33740.7051; %50%//Droop controller parametersfnom = 50; %50.0Vnom = 400; %10.0Pbase = 10000; %500kf1=0.01; %0.01kv1=0.05; %0.05for n=2:500001    %dt = t(n) - t(n-1);  V1d(n)=(2/3)*((sin(theta1(n)))*Va_1(n)+(sin(theta1(n)-((2*pi)/3)))*Vb_1(n)+(sin(theta1(n)+((2*pi)/3)))*Vc_1(n));  V1q(n)=(2/3)*((cos(theta1(n)))*Va_1(n)+(cos(theta1(n)-((2*pi)/3)))*Vb_1(n)+(cos(theta1(n)+((2*pi)/3)))*Vc_1(n));  I1d(n)=(2/3)*((sin(theta1(n)))*Ia_1(n)+(sin(theta1(n)-((2*pi)/3)))*Ib_1(n)+(sin(theta1(n)+((2*pi)/3)))*Ic_1(n));  I1q(n)=(2/3)*((cos(theta1(n)))*Ia_1(n)+(cos(theta1(n)-((2*pi)/3)))*Ib_1(n)+(cos(theta1(n)+((2*pi)/3)))*Ic_1(n)); P_1(n)= V1d(n)*I1d(n) + V1q(n)*I1q(n) ; Q_1(n)= V1q(n)*I1d(n) - V1d(n)*I1q(n) ; %/*  lter */  P_2(n)=(P_1(n)-Pbase);   %//Constant is Ts*wc=50e-6*31.41=0.0015705 0.0015705*  %Pf1old(n)=P_2;  Q_2(n)=(Q_1(n));  % //Constant is Ts*wc=50e-6*31.41=0.0015705 0.0015705*  %Qf1old=Q_2;   %//Droop control m1=kf1*(fnom/Pbase); n1=kv1*(Vnom/Pbase);%theta1(n+1) = (2*pi*Ts)*(f1error-2)   Vd_ref(n) = (Vnom-(Q_2(n)*n1));  %//Note that Vq_ref=0 by default Vq_ref(n)=0; % //PI control takes in Vd_ref and generates Vd_inv and Vq_inv % //PI controller for Vd  error_vd(n)=sqrt(2) * Vd_ref(n) - V1d(n);  prop_vd(n)=kpv * error_vd(n);  int_vd(n)= int_vd(n-1) + kiv * Ts * error_vd(n);  if int_vd(n) > 1      int_vd(n+1) = 1;  elseif int_vd(n) < -1      int_vd(n+1) = -1;  end  id_ref(n)=prop_vd(n) + int_vd(n) - wc * V1q(n);             % //PI controller for Vq  error_vq(n)= -V1q(n);  prop_vq(n)=kpv * error_vq(n);  int_vq(n)=int_vq(n-1) + kiv * Ts * error_vq(n);  if int_vq(n) > 1      int_vq(n+1) = 1;  elseif int_vq(n) < -1      int_vq(n+1) = -1;  end  iq_ref(n)=prop_vq(n) + int_vq(n) + wc * V1d(n);    %//PI controller for Id  error_id(n)=id_ref(n) - I1d(n);  prop_id(n)=kpi * error_id(n);  int_id(n)=int_id(n-1) + kii * Ts * error_id(n);  if int_id(n) > 1      int_id(n+1) = 1;  elseif int_id(n) < -1      int_id(n+1) = -1;  end  Vd_inv(n)=prop_id(n) + int_id(n) - wL * I1q(n);   %//PI controller for Iq  error_iq(n)=iq_ref(n) - I1q(n);  prop_iq(n)=kpi * error_iq(n);  int_iq(n)=int_iq(n-1) + kii* Ts * error_iq(n);  if int_iq(n) > 1      int_iq(n+1) = 1;  elseif int_iq(n) < -1      int_iq(n+1) = -1;  end  Vq_inv(n)=prop_iq(n) + int_iq(n) + wL*I1d(n);   frequency = (fnom-(P_2(n)*m1)); %  0.0003142*         // Constant is 2*pi*Ts=2*3.14*50e-6=0.0003142 %theta1(n+1) = (2*pi*Ts)*(fnom-(P_2(n)*m1)) + theta1(n); %           // Constant is 2*pi*Ts=2*3.14*50e-6=0.0003142 theta1(n+1)=theta1(n)+(2*pi*(frequency*Ts));%if theta1(n+1)>=((2*pi)-(2*pi)*(frequency*Ts))%theta1(n+1)=0;%end %//Inverse dq transform sine_a(n)=Vd_inv(n)*sin(theta1(n))+Vq_inv(n)*cos(theta1(n)); sine_b(n)=Vd_inv(n)*sin(theta1(n))-((2*pi)/3)+Vq_inv(n)*cos(theta1(n)-((2*pi)/3)); sine_c(n)=Vd_inv(n)*sin(theta1(n))+((2*pi)/3)+Vq_inv(n)*cos(theta1(n)+((2*pi)/3));end figure;  subplot(2, 2, 1);  plot(V1d);  title('V1d');  subplot(2, 2, 2);  plot(V1q);  title('V1q');  subplot(2, 2, 3);  plot(I1d);  title('I1d');  subplot(2, 2, 4);  plot(I1q);  title('I1q');   figure;  plot(theta1);  title('theta wt');   figure;  subplot(2, 1, 1);  plot(P_1);  title('Active Power P_1');  subplot(2, 1, 2);  plot(Q_1);  title('Reactive Power Q_1');   figure;  subplot(2, 1, 1);  plot(Vd_inv);  title('Vd_inv');  subplot(2, 1, 2);  plot(Vq_inv);  title('Vq_inv');figure;  subplot(2, 1, 1);  plot(id_ref);  title('id_ref');  subplot(2, 1, 2);  plot(iq_ref);  title('iq_ref');   figure;  plot(sine_a)  title('Inverse DQ')  hold on  plot(sine_b)  hold on  plot(sine_c) figure;   plot(V1d) hold on  plot(V1q)  grid on

Answer 3

[Abdul Moiz Qureshi]

@Sam Chak @Paul Thank you so much for your guidance but here's the twist in this line of code below, that it has to be worked fine with 0 - V1q(n) or - V1q(n) but here it works kinda good when I place 0.089 - V1q(n) instead of zero I still cannot figure it out properly the reason behind it why is there this problem occurs so I'll just continue with this for now as the same condition is applied for all 4 PI here but God knows well what's the matter :).

Once again Thanks to both of you!

error_vq(n)= -V1q(n); 

Answer 4

[Usman]

Do you have the example of multiphase interleaved buck converter close loop?