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

12 vues (au cours des 30 derniers jours)
Abdul Moiz Qureshi le 23 Juil 2023
Réponse apportée : Usman le 30 Juil 2023
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:
• 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));
%//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
##### 4 commentairesAfficher 2 commentaires plus anciensMasquer 2 commentaires plus anciens
Sam Chak le 25 Juil 2023
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 le 25 Juil 2023
Modifié(e) : Abdul Moiz Qureshi le 25 Juil 2023
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 :)

Connectez-vous pour commenter.

### Réponses (4)

Paul le 24 Juil 2023
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?
##### 9 commentairesAfficher 7 commentaires plus anciensMasquer 7 commentaires plus anciens
Abdul Moiz Qureshi le 26 Juil 2023

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

Abdul Moiz Qureshi le 26 Juil 2023
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.

Connectez-vous pour commenter.

Sam Chak le 25 Juil 2023
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).
##### 2 commentairesAfficher AucuneMasquer Aucune
Abdul Moiz Qureshi le 25 Juil 2023
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 le 26 Juil 2023
@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

Connectez-vous pour commenter.

Abdul Moiz Qureshi le 26 Juil 2023
@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);
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

Usman le 30 Juil 2023
Do you have the example of multiphase interleaved buck converter close loop?
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Control and Measurements dans Help Center et File Exchange

R2016a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by